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Note From the Editor 


As previously announced, this is the last issue of The Telecommunications and 
Data Acquisition Progress Report for which printed copies will be produced. Since 
issue 42-121, published on May 15, 1995, The TDA Progress Report has been avail- 
able to all its readers on the World Wide Web at 

http://tda.jpl.nasa.gov / progress_report 

as well as in printed form. Beginning with the next issue, 42-125, due to be pub- 
lished on May 15, 1996, The TDA Progress Report will be published entirely elec- 
tronically at the above-mentioned URL. The TDA Progress Report is linked to 
numerous home pages; among them is the JPL home page at 

http: / /www .jpl.nasa.gov 
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Preface 


This quarterly publication provides archival reports on developments in programs 
managed by JPL’s Telecommunications and Mission Operations Directorate (TMOD), 
which now includes the former Telecommunications and Data Acquisition (TDA) Office. 
In space communications, radio navigation, radio science, and ground-based radio and 
radar astronomy, it reports on activities of the Deep Space Network (DSN) in planning, 
supporting research and technology, implementation, and operations. Also included are 
standards activity at JPL for space data and information systems and reimbursable 
DSN work performed for other space agencies through NASA. The preceding work is 
all performed for NASA’s Office of Space Communications (OSC). 

TMOD also performs work funded by other NASA program offices through and 
with the cooperation of OSC. The first of these is the Orbital Debris Radar Program 
funded by the Office of Space Systems Development. It exists at Goldstone only and 
makes use of the planetary radar capability when the antennas are configured as science 
instruments making direct observations of the planets, their satellites, and asteroids of 
our solar system. The Office of Space Sciences funds the data reduction and science 
analyses of data obtained by the Goldstone Solar System Radar. The antennas at all 
three complexes are also configured for radio astronomy research and, as such, conduct 
experiments funded by the National Science Foundation in the U.S. and other agencies 
at the overseas complexes. These experiments are either in microwave spectroscopy or 
very long baseline interferometry. 

Finally, tasks funded under the JPL Director’s Discretionary Fund and the Caltech 
President’s Fund that involve TMOD are included. 

This and each succeeding issue of The Telecommunications and Data Acquisition 
Progress Report will present material in some, but not necessarily all, of the aforemen- 
tioned programs. 
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The Effect of Aperture Averaging Upon Tropospheric 
Delay Fluctuations Seen With a 
DSN Antenna 

R. Linfield 

Tracking Systems and Applications Section 


The spectrum of tropospheric delay Suctuations expected for a DSN antenna at 
time scales <100 s has been calculated. A new feature included in these calculations 
is the effect of aperture averaging, which causes a reduction in delay fluctuations 
on time scales less than the antenna wind speed crossing time, &5-10 s. On time 
scales less than a few seconds, the Allan deviation o y (A t) oc (At) +1 , rather than 
Oy(At) <x (At) -1 / 6 without aperture averaging. Due to thermal radiometer noise, 
calibration of tropospheric delay Suctuations with water vapor radiometers will not 
be possible on time scales less than fslO s. However, the tropospheric fluctuation 
level will be small enough that radio science measurements with a spacecraft on 
time scales less than a few seconds will be limited by the stability of frequency 
standards and/or other nontropospheric effects. 


I. Introduction 

Radio science experiments with a spacecraft use a one- or two-way radio link between a DSN antenna 
and the spacecraft. The received amplitude and phase are used to measure quantities such as gravitational 
radiation, refractivity profiles of the atmosphere of the planet or planetary satellite near the spacecraft, 
or the structure of a planetary ring system. Any perturbations on the link that are caused by the media 
between the DSN antenna and the spacecraft will corrupt the accuracy of the radio science measurement 
(unless the purpose of the experiment is to study that medium, e.g., the solar plasma). As the link 
frequency used for radio science experiments has increased, the relative magnitude of tropospheric and 
solar plasma phase fluctuations has changed dramatically, because plasma is dispersive at microwave 
frequencies and the troposphere is not. At S-band (2.3 GHz), solar plasma phase fluctuations dominate, 
while at Ka-band (32 GHz), tropospheric phase fluctuations will dominate except at small sun-spacecraft 
angular separations. 

Knowledge of the tropospheric fluctuation spectrum will be useful in planning and analyzing radio 
science experiments. Comparison of fluctuation levels measured by very long baseline interferometry 
(VLBI) and water vapor radiometers (WVRs) has demonstrated that tropospheric delay/phase fluctua- 
tions at microwave frequencies are dominated by fluctuations in water vapor density [l]. 1 The spectrum 


1 C. D. Edwards, “Water Vapor Radiometer Line-of-Sight Calibration Capabilities,” JPL Interoffice Memorandum 
335.1-90-015 (internal document), Jet Propulsion Laboratory, Pasadena, California, March 30, 1990. 
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of tropospheric fluctuations has been measured by VLBI on time scales >20 s [2] and by WVRs on time 
scales >200 s [3]. A model for these fluctuations has been developed [4]. Previous calculations using this 
model (e.g., [5]) have assumed a zero thickness “pencil beam” for a DSN antenna. However, on time 
scales less than «10 s, the nonzero diameter of a DSN antenna will modify this pencil beam spectrum. 
The method for calculating this “aperture averaging” effect is presented in Section II, and results are 
given in Section III. A short summary is given in Section IV. 


li. Calculation of Aperture Averaging 

Both theoretical arguments [6] and observational data [2,4,7] support the idea that tropospheric re- 
fractivity fluctuations obey a Kolmogorov power law. The refractivity structure function D n (r) on scales 
up to at least a few tens of kilometers is 


D n (r ) = ([N( x + r) - N(x )] 2 ^) = C£r 2/3 (1) 

where N(x) is the refractivity at location x (N = n — 1, where n is the index of refraction), and C n 
is the refractivity structure constant. Over time scales up to thousands of seconds, time variations in 
line-of-sight tropospheric delay can be successfully represented by a “frozen flow” model, in which spatial 
variations are convected past the observer by the wind. 

The height dependence of C n is poorly constrained by data. It was modeled by [4] as constant up to 
a 1-km height and zero above that. This “slab height” was subsequently revised to 2 km. 2 Recent WVR 
data from Goldstone [3] show that the observed fluctuations agree with the predictions of this 2-km slab 
model on time scales >400 s but are lower than model predictions on shorter time scales. This result is 
most easily explained by a thickness of the turbulent layer that is significantly larger than 2 km, although 
the data are not adequate to solve for a specific value for this thickness. The effect of a finite thickness 
h of the medium is that fluctuations in the vertical dimension saturate on time scales > h/v w , where v w 
is the wind speed. For time scales At « h/v w , the delay structure function D T oc (At) 5 '' 3 , and for time 
scales At » h/v w , D r oc (A t) 2 ^ 3 [4], The time scales covered in this article are all << h/v w . WVR 
measurements of fluctuations on time scales of 200 s [3] were used to derive the mean C n , once h was 
chosen. For h = 4 km, C n = 3.0 x 10” 8 m -1 / 3 ; these values were used for all calculations presented in 
this article. Over the course of a year, the structure constant at Goldstone exhibits variations of at least 
a factor of 2 about its mean value [3]. 

The tropospheric delay for a pencil beam looking from location x at elevation angle 9 and azimuth 
AZ is 


1 r°° 

Tpencil(x , 6, AZ) = / N [x + r(9, AZ, z)\ dz (2) 

sint/ J o 

where r(9, AZ, z) is the vector from the surface to height z in direction (9, AZ). For a DSN antenna of 
diameter d, the delay is averaged over the circular aperture, to give 

4 r rd/2 

tdsn(x , 9, AZ) = ^2 J J C T pencil [x + F(£, (f>), 9, AZ] d£d(f> (3) 

where r(£, (j>) is the vector in the plane of the aperture, starting from the center, of length £ and at angle 
<f> (the direction chosen for 0 = 0 does not affect the result). Equation (3) assumes that the near- field 


2 R. N. Treuhaft, personal communication, Jet Propulsion Laboratory, Pasadena, California, 1995. 
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DSN beam profile is uniform across its circular cross section. The true profile will be tapered towards 
the edges. As presented in Section III, a change of a factor of two in antenna diameter causes a change 
in the resulting fluctuation level by a factor of «2. Therefore, neglect of beam tapering should cause an 
error of <20 percent in the results reported here. 


A convenient way of calculating and expressing the tropospheric fluctuation spectrum is Allan variance 
<r%(At), defined for a delay process r(t) as [8] 


(At) = 


^[r(f + 2A t) — 2r(f + At) + r(£)] 2 ^ 
2 (At) 2 


( 4 ) 


Expanding Eq. (4) and assuming that averaged quantities are independent of time and position gives, for 
the tropospheric delay fluctuations measured by a DSN antenna, 


_ 2 / a ja 3 ( t dsjvW) 4 (T DSN (t + At)r DSN (t)) , (TDSN(t + 2At)r DSN (t)) , rS 

a y {M) -~(AW (At)* + (At)* (5j 

We can evaluate this expression with Eq. (3) and the frozen-flow assumption. Taking the middle term in 
Eq. (5) as an example, 

jg p2it p2ir pd/2 pd/2 

(T D SN(t + At)T DSN (t)) I J J o J q (r P encii[x + v w At + r(Z, </>)} T pencil [£ + r(? , </>')]) 

x d£d£' d(j)d<p' 


It has been assumed that averaged quantities are position independent, so that 


(^"pencil (^ 1 7 AZ)Tpe nc n(X2' t AZ)) (Tp enc n(x , 6 , A.Z)*) ^ AZ) Tp enc n(x 2 , 0 j A.Z)} ^ 

= ( T Zencil(x,0,AZ))-^D r ( |fi-X 2 |) 


where D r is the pencil-beam delay structure function, defined as 

Dt (? ) = \Tpcnctl i x AZ) Tpencil(%i@ ■> AZ ) ] ^ 

In evaluating Eq. (5), the ( Jp enc u (^ terms add to zero, and 0 y(At) can be expressed as a four- 
dimensional integral of three D T terms: 
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•ig p 2 tt p2i r pd/2 /*d/2 

wi M l «J. { n « 


^=R«»-nr,«i 


B - \v w A t + r(£, 0) - r(£', <£')| 


C = \2v w At + r(f , <f>) - f(£', ^') 


-^D T (A) + 2D T (B)-iD T (C) 


> ( 6 ) 


The horizontal distance scales are d and v w A t. As a result, cr y (A t) will be suppressed, relative to its 
d = 0 value, for values of At smaller than a few times the wind speed crossing time of t croas = d/v w . 
For a typical wind speed in the lower few kilometers of 8 m/s and for d — 34 m, t croas « 5 s. For 
At « t crosa , ay(At) will be much lower than its d = 0 value, due to the averaging of many small-scale 
delay fluctuations over the much larger DSN aperture. 

For any DSN tracking scenario at an observing frequency high enough (>8 GHz) for tropospheric 
delay fluctuations to dominate over those from plasma, the Earth’s troposphere will lie entirely within 
the near field of the antenna. The effective tropospheric volume sampled by the DSN antenna will then 
be cylindrical. For tracking of planetary spacecraft, the antenna will move at a nearly sidereal rate. This 
tracking will cause an effective velocity v e f / through the troposphere, at height 2 and elevation angle 6, 
of 


Veff = 


Clz 

sin# 


The Earth’s rotation rate is f2. For z = 2 km (half the turbulent slab height), v e ff = 0.15/ sin 9 m/s. 
T his is much less than the wind speed of 8 m/s used in the calculations. Sidereal tracking, therefore, will 
not significantly modify the results presented in this article. 

The Fresnel length scale l for observations at a wavelength A with phase perturbations occurring at a 
distance L is l « V\L. For A = 1 cm and L — 2 km (the middle of the turbulent troposphere), / « 4 m. 
On smaller scales, geometric optics is not completely valid. The results for time scales <4 m/s/v w = 0.5 s 
should, therefore, be regarded as approximate. 


III. Results 

Equation (6) was evaluated numerically, for a range of time intervals, and for both 34-m- and 70-m- 
diameter antennas. A wind speed of 8 m/s was assumed. The results for antennas pointed in the zenith 
direction are shown in Fig. 1. Values of cr y (At) for a pencil beam antenna are shown for comparison. 

For time intervals shorter than 10 s, the Allan deviation for measurements made with a DSN antenna 
are less than those made with a pencil beam. For time intervals less than «2 s, cr y (At) oc (A t) k , where 
k « +1. This proportionality can be roughly understood as follows: The effect of fluctuations over a time 
interval At are confined to a region of size wu w At. For v w A t « d, there are N regiona ~ (d/(v w At)) 2 
such regions over the antenna aperture. Averaging the delay fluctuations of these regions will cause a 
reduction by a factor of ~ N regions ~ d/{v w At). The slope of the pencil beam a y (At) is w— 1/6, so the 
slope of the DSN cr y (At) will be w+5/6. 

As At decreases, the amount of computer time needed to calculate a y (At) increases rapidly. For small 
At, the terms A, B , and C in Eq. (6) are nearly equal, and high numerical precision is needed in order 
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Fig. 1. The Allan deviation calculated for a pencil 
beam antenna and for antenna diameters of 34 and 
70 m is shown as a function of time interval for 
typical Goldstone conditions and the zenith 
direction. The Allan deviation at other elevation 
angles will be somewhat larger, increasing by a 
factor of 1.8 at a 10-deg elevation angle. The thermal 
noise from a WVR with the characteristics expected 
for the Cassini radio science calibration system is 
included. 


to determine the differences in the D r terms. Furthermore, all derivatives of D n and all but the first 
derivative of D r become singular at zero separation, so that numerical integration methods have difficulty. 


Calculations for At as small as 0.02 s were possible in the zenith direction, because D T (r) could be 
evaluated separately, with an empirically fitted analytical formula used for D T in Eq. (6) to reduce the 
number of dimensions in the integration from 6 to 4. Calculations at other elevation angles were limited 
to a narrower range of time intervals (1-100 s). These calculations showed that cr y (At) increases slowly 
with decreasing elevation angle. Relative to its value at the zenith, a y (At) was 1.16 times larger at a 
30-deg elevation angle and 1.82 times larger at a 10-deg elevation angle. 


All cr y {At) values in Fig. 1 are linearly proportional to the structure constant, C n . Therefore, the 
actual cr y (At) values at a DSN site will range up and down from those in Fig. 1 by a factor of at least 2, 
depending upon the season, time of day, and weather conditions. Winter nights tend to have the lowest 
C n values, with summer days having the largest values. 


Tropospheric delay fluctuations at microwave frequencies can be calibrated by WVRs, which measure 
thermal emission from water vapor in the vicinity of its 22-GHz spectral line [9]. However, on short 
time scales, the thermal noise from a WVR becomes the limiting error source. The main vapor-sensing 
channel of current WVRs is at a frequency near 21.5 or 23.8 GHz, where a 1-K brightness temperature 
corresponds to «6 mm /c = 2 x 10“ 11 s of path delay (c is the speed of light). The thermal path delay 
equivalent noise in a WVR measurement of integration time for a total system temperature T sys and 
bandwidth BW is 


N(t int ) = 


•-sys 


2 x 10- 11 s 2 x 10 -13 Ti 


100 


VBWtTr, 


IK 


BW™ 


(7) 


where Tioo = T ay3 / 100 K, BW wo = BW / 100 MHz, and t int3 = t int /l s. 
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To set a firm lower limit on the effect of WVR thermal noise, we assume that At = t int (i.e., we give 
up all information on time scales shorter than At). Future uncooled WVRs currently under design are 
expected to have T sys « 300 K and BW « 400 MHz. 3 The Allan deviation from thermal WVR noise for 
A t 3 = At/1 s (note that a y (At) = iV(f; nt )\/3/Af for thermal noise) then is 


CTy(At) 


5 x HT 13 

A tp 


(WVR thermal noise) 


( 8 ) 


This WVR thermal noise curve is plotted in Fig. 1. If we require a y (At) from WVR thermal noise to be 
at least 3 times smaller than the tropospheric cr y (At) in order to achieve useful calibration, then, for any 
elevation angle >10 deg, the minimum time interval A t m i n for WVR calibration is 


At 


min 


10 s 


Because the slope of the tropospheric cr y (At) is small at At ~ 10 s and the slope of the WVR thermal 
noise <J y (A t) is steep, the minimum useful calibration interval is relatively insensitive to the values of 
T sys and BW. On time scales less than «10 s, tropospheric fluctuations cannot be calibrated with any 
known technique. 


IV. Summary 

The effect of aperture averaging upon radio science measurements is that tropospheric fluctuations will 
not be important on short time scales. For occultation measurements of Saturn’s rings with Cassini, the 
shortest time scale of interest is the desired physical resolution («100 m) divided by the spacecraft orbital 
velocity («10 km/s), or ~0.01 s. 4 These measurements will use an onboard oscillator with a y (At) > 
10 -13 . Even if future flight oscillators were a factor of 10 more stable than the one on Cassini, the 
troposphere would be more stable than the oscillator for At < 1 s. The tropospheric fluctuations on 
At < 10 s appear to form a fundamental lower limit for radio science measurements using an Earth 
antenna, because WVR thermal noise precludes calibration of the fluctuations. 
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An Overview of the GOLD Experiment Between 
the ETS-VI Satellite and the Table 
Mountain Facility 

K. E. Wilson 

Communications Systems and Research Section 


The Ground /Orbiter Lasercomm Demonstration is a demonstration of optical 
communications between the Japanese Engineering Test Satellite (ETS-VI) and an 
optical ground transmitting and receiving station at the Table Mountain Facility 
in Wrightwood, California. Laser transmissions to the satellite are performed for 
approximately 4 hours every third night when the satellite is at apogee above Table 
Mountain. The experiment requires the coordination of resources at the Communi- 
cations Research Laboratory (CRL), JPL, the National Aeronautics and Space De- 
velopment Agency (NASDA) Tsukuba tracking station, and NASA’s Deep Space 
Network at Goldstone, California, to generate and transmit real-time commands 
and receive telemetry from the ETS-VI. Transmissions to the ETS-VI began in 
November 1995 and are scheduled to last into the middle of January 1996, when 
the satellite is expected to be eclipsed by the Earth’s shadow for a major part of 
its orbit. The eclipse is expected to last for about 2 months, and during this period 
there will be limited electrical power available on board the satellite. NASDA plans 
to restrict experiments with the ETS-VI during this period, and no laser transmis- 
sions are planned. Posteclipse experiments are currently being negotiated. GOLD 
is a joint NASA-CRL experiment that is being conducted by JPL in coordination 
with CRL and NASDA. 

L Introduction 

The Ground/Orbiter Lasercomm Demonstration (GOLD) is a joint NASA/Communications Research 
Laboratory (CRL) optical communications experiment to evaluate one-way and two-way optical commu- 
nications under a range of atmospheric conditions and to demonstrate optical ranging. GOLD’S objectives 
and goals are outlined below. 

The objectives are to 
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(1) Demonstrate two-way spatial acquisition /tracking of laser beams with a spacecraft 

(2) Accomplish one-way and two-way optical data transfer to a spacecraft and measure bit- 
error rates 

(3) Accumulate 10 elapsed hours of transmission/reception experience and 30 Gbytes over a 
6-month period 

(4) Compare downlink atmospheric transmission losses with similar data from the Table 
Mountain Facility’s (TMF’s) atmospheric visibility monitoring (AVM) observatory 

The goals are to 


(1) Gather atmospheric transmission statistics from an actual space- to-ground link and com- 
pare them with corresponding statistics from the AVM system 

(2) Build a database for optical link acquisition/reacquisition times 

(3) Validate optical communications-link performance prediction tools 

(4) Demonstrate optical ranging to 10-m accuracy 

GOLD experiments will use the 0.6-m and 1.2-m telescopes located at NASA’s TMF to communicate 
with the Japanese Engineering Test Satellite (ETS-VT). The experiment concept is depicted in Fig. 1. An 
argon-ion laser coupled to the 0.6-m telescope transmits a 1.024-Mbps Manchester-coded pseudorandom 
noise (PN) sequence to the spacecraft. The ETS-VI in turn uses its GaAs laser to transmit a similar PN 
sequence to the 1.2-m ground receiver located approximately 60-m from the transmitter site. 


830-nm DOWNLINK 


1.2-m RECEIVER 



51 4.5-nm UPLINK 


0.6-m TRANSMITTER 


Fig. 1. Conceptual drawing of GOLD. The transmitter at TMF uplinks a 51 4.5-nm com- 
munications signal to the ETS-VI. The satellite downlinks a 830-nm 1-Mbps signal. 
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The ETS-VT was launched into orbit on August 28, 1994. It was originally intended to be in a 
geostationary orbit above Japan, but difficulty with one of its motors has resulted in the satellite now 
being in a geotransfer orbit. To make maximum use of the spacecraft’s subsystems in its current orbit, 
researchers at CRL have encouraged both NASA and European Space Agency (ESA) experimenters to use 
the optical communications subsystem on board the ETS-VI. In response to this, the National Aeronautics 
and Space Development Agency (NASDA) has refined the satellite’s orbit to facilitate the use of the laser 
communications equipment (LCE) by experimenters at JPL. Yet, because the satellite’s elliptical orbit 
takes the satellite through the Van Allen belts, its power generation capabilities have been degraded 
and its expected life reduced. Between mid-January and mid-March 1996, the satellite will go into a 
2-month-long eclipse. During this time, very little power will be generated on the spacecraft, and only low- 
power experiments will be feasible. It is unclear whether optical communications falls into this category 
of experiments. Posteclipse experiments are not scheduled as yet because of the large uncertainty in 
whether the satellite batteries will survive the long eclipse. 


This article describes the work on the GOLD project to date. The three phases of the GOLD ex- 
periment are described in Section II, along with a brief description of the satellite predict generation. 
The data link between CRL and TMF that is used to confirm the acquisition of the uplink beam and to 
communicate with CRL experimenters is also described in Section II. The transmitter and receiver optics 
and electronics are described in Sections III and IV, respectively, and the AVM station that measures the 
atmospheric transmission at optical wavelengths is described briefly in Section V. Samples of the data 
collected to date are presented and discussed in Section VI. Conclusions are in Section VII. 


II. GOLD Preparation and Operations Phase 

A. GOLD Operations and Experiment Phases 

The ETS-VI satellite is in a recurrent orbit that makes it visible every 3 days at night from TMF. The 
GOLD schedule given in Table 1 extends over the period from October 30, 1995, to January 13, 1996. 
The table shows the laser communications and associated Deep Space Network (DSN) satellite control 
times. The DSN begins satellite control 2 hours before the laser transmission times and ends 90 minutes 
after laser transmission. 


Table 1 shows the transmission times of the GOLD experiments. The experiment is broken into 
three phases over the period from November 1995 to January 1996. Phase 1 is a beacon uplink with a 
1.024-Mbps downlink detection. The goals of this phase are to 


(1) Complete integrated system testing 

(2) Measure the uplink optical beam divergence at the satellite 

(3) Evaluate the benefits of spatial diversity of the uplink beacon, i.e., single- and dual-beam 
beacon transmissions 

(4) Determine the appropriate uplink power levels to maintain spacecraft closed-loop track- 
ing of the ground station 

(5) Evaluate the spacecraft’s ability to point to the TMF ground station, and develop ap- 
propriate ground-station acquisition strategies 

(6) Establish operations procedures for the phases discussed below 
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Table 1. Laser transmission times for the GOLD experiment. 


Date 

Start, 
time UT 

End, 
time UT 

Duration, 

h:min 

30 October 1995 

11:55 

15:00 

2:05 

2 November 1995 

11:45 

15:00 

3:15 

5 November 1995 

11:30 

14:50 

3:20 

8 November 1995 

11:20 

14:40 

3:20 

11 November 1995 

11:15 

14:30 

3:15 

14 November 1995 

11:05 

14:10 

3:05 

17 November 1995 

11:05 

14:10 

3:05 

20 November 1995 

10:45 

14:00 

3:15 

26 November 1995 

10:35 

13:30 

2:55 

29 November 1995 

9:52 

13:55 

4:01 

2 December 1995 

9:36 

13:50 

4:14 

5 December 1995 

9:20 

13:47 

4:27 

8 December 1995 

9:04 

13:43 

4:39 

11 December 1995 

8:48 

13:40 

4:52 

14 December 1995 

8:35 

13:35 

5:00 

17 December 1995 

8:25 

13:24 

4:59 

20 December 1995 

8:14 

13:13 

4:59 

23 December 1995 

8:04 

13:03 

4:59 

29 December 1995 

7:44 

12:41 

4:57 

4 January 1996 

7:24 

12:20 

4:56 

7 January 1996 

7:14 

12:09 

4:56 

10 January 1996 

7:04 

11:59 

4:55 

13 January 1996 

6:54 

11:48 

4:54 


Phase 2 extends over the month of December 1995, and the goals are to demonstrate two-way optical 
communications, measure the bit-error rate (BER) on the uplink, and evaluate the spacecraft’s ability to 
regenerate an uplinked data sequence. 

Phase 3 extends from January 4, 1995, to January 13, 1996. The primary goal during this period is 
to use the spacecraft’s uplink signal regeneration capabilities to measure the one-way light-time to the 
spacecraft. These data will be processed through JPL’s orbit determination program (ODP) to generate 
an element set to predict the spacecraft ephemeris. 

Table 2 gives a brief description of the results of each night’s operations as of the time of writing. Once 
the required LCE gimbal angle relative to the satellite’s orientation was established on November 11, 1995, 
the uplink was acquired within seconds of transmission on subsequent nights. 

B. Satellite Predicts Generation 

The satellite predicts for the TMF telescopes are generated from a satellite osculating-element set that 
is faxed to JPL from NASDA’s Tsukuba tracking facility. Ephemeris files are generated at JPL for the 
transmitter and receiver telescopes, separately with offsets to account for the point ahead angle. The 
files are electronically transmitted to TMF and are loaded into the telescope control program (TCP). 
The TCP uses a spline-fitting routine to interpolate between the ephemerides to generate a “smooth” 
telescope-pointing file to track the ETS-VI satellite. 
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Table 2. Brief description of results for each night’s operations. 


Date 


Description 


30 October 1995 
2 November 1995 
5 November 1995 
8 November 1995 


11 November 1995 


14 November 1995 
17 November 1995 


20 November 1995 
23 November 1995 
26 November 1995 


Heavy cloud cover precluded transmission 
Transmitted laser beacon to satellite; no uplink detected 
Transmitted laser beacon to satellite; no uplink detected 
First detection of uplink beacon 
Spacecraft gimball was being scanned 
Downlink detected and recorded 

Uplink detected for 45 min; became sporatic thereafter 
Spacecraft has difficulty tracking the ground station 
Downlink detected and recorded 
Uplink detected, and downlink data recorded 
Spacecraft detector gains adjusted to prevent saturation 
Satellite acquired within seconds of initiating beacon 
Uplink detected for 3 h 
0.6-m telescope malfunction 
Thanksgiving holiday 

Heavy cloud cover over TMF; no transmission 


C. Data Links 

GOLD requires the simultaneous real-time coordination of experimenters at NASDA, CRL, and the 
DSN to accomplish the laser transmission from TMF. The command and data flow diagram is shown 
in Fig. 2. In addition, during the early part of the experiment, the Federal Aviation Administration’s 
(FAA’s) Air Traffic Control Center in Palmdale, California, is assisting in aircraft avoidance strategies 
until an aircraft detection radar system at TMF is installed. 


Commands to control the satellite’s attitude during the experiment are generated at NASDA. CRL 
generates commands to control the LCE and sends them to NASDA. These commands are then trans- 
mitted from NASDA to DSS 27 at Goldstone, which relays them to the satellite. 


The data link between CRL, JPL, and TMF provides near-real-time feedback on the measured onboard 
laser communications equipment sensors to the experimenters at TMF. The delay between the spacecraft 
transmission and the TMF reception is approximately 15 s. Spacecraft attitude control system (ACS) 
data are transmitted along with LCE data via S-band (2.3 GHz) telemetry to DSS 27. The DSN transmits 
these data to the NASDA Space Center at Tsukuba, where they are demodulated. The LCE data are 
transmitted to CRL for processing. CRL in turn transmits an ASCII data stream showing time, charge- 
coupled device (CCD) level, quadrant detector (QD) level, spacecraft laser-diode bias current, and coarse 
and fine tracking-sensor errors via integrated digital network services (ISDN) link to JPL. Because there 
are no ISDN links to Wrightwood, a 56-kbps telephone link is used to transmit the data from JPL to 
TMF. 


IIS. The Transmitter 

The transmitter consists of an argon-ion laser coupled to the 0.6-m telescope at TMF. The telescope 
is located in Building TM-12 at TMF, and its surveyed position is as follows: 
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Parameter 

Position 

Longitude 

117° 40' 52.55" 

Latitude 

34° 22' 53.49" 

Altitude 

2.286 km 


The telescope is used in the coude mode, which allows light to be coupled from large, high-power lasers 
into the telescope. The uplink laser is a prism-tuned coherent Innova-100 argon-ion laser that delivers 
a maximum 14.5-W linearly polarizer laser-light output power at 514.5 nm. At maximum power, the 
laser output is multimode. The laser is typically operated at an output power of about 13 W to achieve 
good beam quality. Coalignment of the laser beam with the telescope axis was achieved by adjusting 
the position of the telescope’s secondary mirror until the focus was brought to a position on the optical 
bench. At this setting, the telescope’s focal ratio was f/41, i.e., a focal length of 26.4 m. 

A schematic of the optical train is shown in Fig. 3. A Conoptics electro-optic modulator impresses 
the uplink data stream on the optical carrier. The modulator consists of four potassium-dihydrogen- 
phosphate crystals and a polarizer. A data formatter — a Firebird 6000 bit-error rate tester (BERT) is 
used to generate a basic data pattern that is amplified — generates 0- to 1-V square-wave modulation 
that is amplified to the modulator’s half-wave voltage in the driver, with 0 V corresponding to maximum 
transmission through the modulator. 

After modulation, the beam is incident on a concave/convex lens pair that sets the beam divergence 
out of the telescope. This is nominally set at 20 prad. A beam splitter separates the beam into two 
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Fig. 2. Real-time coordination of several organizations is needed to accomplish the GOLD experiment. NASDA 
transfer commands to the DSN, which uplinks them to the ETS-VI. The spacecraft attitude and LCE sensor status 
are downlinked via S-band telemetry to Goldstone. The telemetry is demodulated by NASDA and forwarded to CRL 
for processing. During the optical uplink, CRL transmits LCE sensor data to JPL and TMF. 
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Fig. 3. Schematic of the optical train for the GOLD experiment 
showing the laser, the modulator, and the optical beam delay system 
used to provide temporal and spatial diversity of the optical beams 
transmitted to the satellite. Also shown are the detector used to 
monitor the modulated wave form and the CCD camera used to 
image the satellite. 


equal parts, one of which goes through a 25-cm optical delay line with a path length difference greater 
than the laser’s coherence length; A 2 /AA is 10 cm. Both beams are reflected from a high-power dichroic 
beam splitter and are brought to a focus at the iris that is located at the f/41 focus of the telescope. 
From there the beams diverge and are reflected by the third coude flat and into the telescope. The beams 
are made incident on opposite sides of the 0.6-m telescope’s primary mirror, a distance greater than the 
size of an atmospheric coherence cell. The use of spatial and temporal diversity mitigates the effects of 
atmospheric scintillation on the uplink beacon, thereby allowing the satellite’s tracking system to better 
track the uplink beacon. 

The CCD was a Pulnix camera and image intensifier with sensitivity down to 10 -6 lux. This high- 
sensitivity camera has enabled us to track the satellite around apogee, where its brightness has varied 
from that of a magnitude 12 to a magnitude 14 star, depending on the phase angle of the solar panels. The 
avalanche photodiode (APD) transmitter detector monitors the modulation of the transmitted signal. In 
the ranging phase of the GOLD experiment, this detector will monitor the modulation sequence sent to 
the satellite, and these data will be transmitted to the receiver facility to initiate the timing sequence. 


IV. The Receiver Facility 

The optical receiver weighs approximately 30 kg and is mounted at the flange of the 1.2-m (f/29.5) 
telescope’s bent Cassegrain focus. It consists of two CCD cameras and a 3-mm diameter low-noise APD 
(see Fig. 4). The cameras are a wide-field Cohu with an image intensifier for satellite acquisition and 
tracking and a Spectra Source CCD for atmospheric seeing measurements. These detectors are coaligned 
on an optical bench assembly to ensure that the downlink transmission is incident on both the tracking 
and communications detectors. 

Satellite acquisition at the receiver is accomplished in a series of steps. It begins with calibration of the 
telescope’s pointing direction relative to a 0.25-m Meade and a 0.4-m guiding telescope attached to the 
telescope’s frame. A bright calibration star is acquired in the wider-field Meade and is transferred to the 
satellite CCD tracking detector in the smaller 1.2-m telescope field. Because the blind point accuracy of 
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Fig. 4. Schematic of the optical receiver located at the focus of the 1.2-m 
telescope. CCD detectors in the optical train track the satellite and measure 
atmospheric seeing. The APD detects the 1.024-Mbps optical downlink data 
stream. 

the telescope depends on separation between the calibration star and the satellite positions, the telescope 
is moved to a star in the vicinity of the satellite and the pointing offsets are noted for calibration. Near 
apogee, the satellite brightness ranges from the 12th to the 14th magnitude. Because of the large dynamic 
range between the third- to fifth-magnitude calibration stars and the satellite, an optical density no.-2 
filter is placed in an electronically switchable holder located in front of the tracking camera. The filter is 
placed in the optical beam during the calibration and is switched out to acquire the satellite. 

Atmospheric seeing measurements provide essential data for evaluating the optical link performance. 
The predictions of theoretical models that incorporate scintillation effects into the link performance will 
be compared with experimental results to validate the models. The data are taken at 15-min intervals 
and are reduced and transferred to the data analyst at the end of each night’s run. 

Downlink data-recovery electronics starts with the APD amplifier and signal conditioner. Depending 
on the mode of operation, the downlinked data are recovered through a bit synchronizer to recover the 
clock and then stored on a digital data recorder. This is shown in Fig. 5. When operating in the ranging 
or regeneration modes, before storage, the data are correlated with the transmitted data from the BERT 
located in the transmitter facility. 


V. AVM Observatory 

The atmospheric attenuation during the experiment is measured using the AVM observatory located 
at TMF. The AVM measures and records the intensity of selected stars over the course of the experiment. 
Star intensities are measured using three Johnson standard [V (100-nm wide centered at 560 nm), R 
(200-nm wide centered at 700 nm), and I (200-nm wide centered at 860 nm)] astronomical filters and 
three 10-nm-wide interference filters centered at the 532-nm, 860-nm, and 1064.2-nm wavelengths. The 
attenuation measurements made at 532 nm and 860 nm are used to estimate the atmospheric trans- 
mission at the uplink (514.5-nm) and downlink (830-nm) laser wavelengths. Transmission data at these 
wavelengths are calibrated using data taken through the astronomical V and I filters. 


VI. Data Recovery and Data Processing 

Data recovery from the satellite is accomplished using both S-band and optical downlinks. The up- 
link signal power is sampled at the CCD and QD at 1 Hz and transmitted on S-band telemetry. A 1-h 
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Fig. 5. Schematic of downlink data recovery. The APD output is bit synchronized and recorded for 
later processing. The 602-A digital oscilloscope monitors the low-pass-filtered downlink signal for 
long-term signal amplitude variations. 
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Fig. 6. A 1-hour sample of the uplink beacon as measured by the onboard CCD and QDs: (a) satellite 
coarse-tracking CCD sensor and (b) satellite fine-tracking QD sensor. The data show that the uplink 
signal is strong enough to saturate these detectors, in subsequent experiments, improved tracking 
performance was achieved by reducing both the transmitted intensity and the detector gain. 

section of the uplink laser signal measured on the November-14 pass is shown in Fig. 6. The data clearly 
show saturation of the uplink power measurement on both the CCD and the QD. The coarse-pointing 
tracking system that is servoed around the CCD camera acquires the uplink first. The detected signal 
is then centered in the camera’s field of view and tracked by the fine-tracking loop servoed around the 
QD’s output. The delay between acquisition and tracking is on the order of minutes and can be clearly 
discerned in the two traces presented in Fig. 6. The dropouts seen in the data file are due to a combination 
of scintillation on the uplink beam and satellite and LCE pointing. 

The LCE can impress three different data types on the downlink optical carrier. The first is the PN 
mode, in which an onboard PN sequence is transmitted to the ground station. This is shown in Fig. 7, 
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where the random bit flips of the modulation that make up the 1.024-Mbps PN data stream are clearly 
evident. The second is the telemetry mode, where the LCE data are transmitted at 128 kbps (with 
redundant transmission and Manchester modulation up to the 1.024-Mbps channel rate). Included in 
this data stream are the laser diode levels, the CCD and QD signal levels, the fine- and coarse-tracking 
levels, the APD communication detector signal levels, and the measured uplink BER data. A sample of 
this data stream is shown in Fig. 8. The data show bit flips occurring in multiples of eight, consistent 
with the x8 multiplication of the 128-kbps data stream to achieve the 1.024-Mbps data rate. 

The third mode, the regeneration mode, is shown in Fig. 9. In this mode, a 1-MHz square-wave uplink 
is detected by the communications detector and then retransmitted back to the ground station. Figure 9 
clearly shows the regenerated 1-MHz uplink. It also shows the signal fades at a frequency inconsistent 
with atmospheric effects. We believe that these effects are caused by distortion in the modulated uplink 
beam. The causes of this distortion are being investigated. This regeneration feature will be used to 
perform the turnaround optical ranging experiments. 


VII. Conclusions 

GOLD is an international cooperative experiment that has demonstrated two-way optical communi- 
cations between a satellite at a geostationary distance and an optical ground receiver. GOLD’S thrust is 
to measure and understand the performance of the two-way optical link under a variety of atmospheric 
attenuation and turbulence conditions. The data accumulated from this experiment will enable better 
definition of the performance of optical communications systems for mission designers. Satellite rang- 
ing by detection and retransmission of an uplinked code is a time-tested approach in rf communications 
systems. A key goal of GOLD is to demonstrate this capability at optical frequencies for the first time. 

To date we have amassed several gigabytes of data and videotapes of the optical links, along with 
measurements on the uplink scintillation. These data will be processed and the findings reported in 
future articles. 



Fig. 7. Sample of 1.024-Mbps Manchester-coded PN sequence downlink 
telemetry from the LCE showing random bit flips. 
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Fig. 8. Downlinked satellite telemetry at 1 28 kbps. The data show bit flips in 
multiples of eight consistent with 8X repetition of the bit pattern. 
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Uplink data from recent free-spaee optical communication experiments carried 
out between the Table Mountain Facility and the Japanese Engineering Test Satel- 
lite are used to study fluctuations caused by beam propagation through the atmo- 
sphere. The influence of atmospheric scintillation, beam wander and jitter, and mul- 
tiple uplink beams on the statistics of power received by the satellite is analyzed and 
compared to experimental data. Preliminary analysis indicates the received signal 
obeys an approximate lognormal distribution, as predicted by the weak-turbulence 
model, but further characterization of other sources of fluctuations is necessary for 
accurate link predictions. 


I. Introduction 

An optical communication link between an Earth-orbiting satellite and a ground-based receiver involves 
propagation of a laser beam through the atmosphere. Spatial and temporal variations in the index of 
refraction of the air that forms the atmosphere severely degrade the quality of the beam. These variations, 
typically lasting on the order of one thousandth of a second at optical and near-infrared (NIR) frequencies, 
in the index of refraction result in random changes in the amplitude and phase of the arriving wave. The 
amplitude and phase fluctuations manifest themselves in effects such as scintillation, beam broadening, 
and beam motion or wander. Scintillation can be thought of as interference between partial waves 
propagating through different paths (or turbulent cells), which results in fades (destructive interference) 
or surges (constructive interference) at the receiver. Deep fades causing signal loss or strong surges causing 
saturation of the quad-detector can force loss of track. On the uplink, with the exception of heavy cloud 
cover, atmospheric scintillation is usually the limiting factor in an optical link. For downlink, however, 
the large size of the receiver’s aperture typically compensates for scintillation by averaging over fades and 
surges. Beam motion, primarily on the uplink, can often be considered as atmospheric-induced jitter in 
the pointing of the laser beam. Such motion of narrow beams can cause deep fades in the signal due to 
the Gaussian nature of the spatial beam profile. 

In this article, we use data from the recent Ground/Orbiter Lasercomm Demonstration (GOLD) 
experiments to study the atmosphere-induced fluctuations in signal power on the uplink. The experiment’s 
objective is to establish a 1-Mbps optical communication link between the laser communication experiment 
(LCE) package on board the Japanese Engineering Test Satellite (ETS-VI) and the Table Mountain 
Facility (TMF) near the Jet Propulsion Laboratory (JPL). The link is established by first transmitting a 



beacon from TMF to ETS-VT using a priori knowledge of the satellite’s orbit. The beacon is subsequently 
acquired and tracked by the LCE package on the satellite. Once the tracking loop is activated on the 
LCE, an onboard laser is turned on for downlink. The downlink beam is then picked up by the receiver 
at TMF to complete the link. Once the link is established, both the downlink laser and uplink beacon 
can be modulated for transmission of data. Here, we shall restrict our analysis to scintillation-induced 
signal fluctuations on uplink and ignore the communication aspect of the experiment. 

In the following analysis, we compare the statistics of the power received by the detectors on board the 
satellite when a continuous wave (CW) beacon is transmitted from TMF to the probability distribution 
function (pdf) predicted by existing models for atmospheric propagation. In Section II, we begin with a 
brief summary of the theory that takes into account atmospheric scintillation, jitter, and multiple uplink 
beams. Section III provides the specifics of the experiments and a typical link budget for the uplink. 
The experimental data are presented in Section IV. The final section discusses the observed results in 
comparison with theoretical predictions and offers potential sources of discrepancy between the theory 
and the experiment. 


II. Theory 

Consider a communication channel between a ground station and an orbiting satellite that is specific 
to GOLD, as shown in Fig. 1. A detailed discussion of the influence of the atmosphere on the statistics 
of the channel during beacon uplink follows. Figure 2 shows a system-level representation of the optical 
channel. We will show that the power, Pr, received by the photodetector on the satellite is the product 
of three quantities: 


Pr = PoIS 


where Po is the power received in the absence of a turbulent atmosphere, I is a random variable with 
a beta-distribution caused by atmospheric and pointing jitter, and S typically is a lognormal random 
variable due to scintillation. The maximum receivable power, Po, at the satellite is assumed to have zero 
variance. Thus, fluctuations due to the laser and detector noise are assumed to be negligible compared to 
scintillation and beam jitter. In the following sections, the major contributions to the observed statistics 
of the received power are described. 
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Fig. 1. Geometry of a ground-to-orbiter laser 
communication link showing transmission 
from TMF to ETS-VI. 
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Fig. 2. Components and channels in an optical communication system. 

A. Free-Space Propagation 

For an optical transmitter with continuous wave (CW) laser power of Pr W, the power received by 
the detector on board the satellite in the absence of atmospheric turbulence is primarily governed by 
free-space-propagation loss and optical absorption. Taking into account optical losses at the transmitter 
and receiver and atmospheric absorption, we find that the received power is given by 

D sec0 ^ ArPt 

Po — VroeVtoeVo . _ .2 

* ( ze bd ) 

where Z is the communication range in meters; is the laser beam divergence in radians; 6 is the 
zenith angle to the receiver looking from the transmitter station; % oe is the transmitter optical efficiency; 
T] roe is the receiver optical efficiency; r/ 0 is the atmospheric transmission at zenith; and Ar is the area of 
the receiver’s aperture in square meters. The exponent sec 9 is often referred to as the air mass. We may 
to a good approximation assume that all quantities in the above equation are constant over the duration 
of the experiment, which typically lasts a few hours. 

B. Light Wave Propagation in Turbulent and Random Media 

Wave propagation through lossless random media, like the atmosphere, is often characterized by the 
structure constant C%, which is a measure of the variance of the index of refraction, n, due to inhomo- 
geneities in the medium [1] . For modeling the atmosphere, one often needs to know the altitude variation 
of the structure constant and its dependence on weather conditions. Though modeling the atmosphere 
is a difficult task, there exist several empirical and parametric models for the altitude variation in the 
structure constant of the atmosphere. The most common of these is the theoretically based Hufnagel 
model, in which the dominant turbulence arises from winds in the 5- to 20-km altitude range [2]. The 
altitude of TMF (2.2 km above sea level) is well outside the range of validity of the Hufnagel model 
(over 5-km above sea level). The validity of the improved Hufnagel-Valley model [3], on the other hand, 
extends down to the ground. This model, however, predicts a larger than observed turbulence at TMF 
by overestimating the low-altitude contribution. We, therefore, use the experimentally based Air Force 
Geophysical Laboratory ( AFGL) CLEAR I night model [4] . The AFGL CLEAR I night model is based on 
the average of a large number of measurements of the atmospheric properties in the New Mexico Desert. 
The model provides the functional dependence of C\ on altitude down to 1.23 km above mean sea level: 

log 10 C„ = A + Bh + Ch 2 + D exp 



where h is the altitude in km above mean sea level. The constants A , B, C, D, E, and F are defined for 
three different ranges of altitude and listed in Chapter 2 of [3] . Most quantities of interest are usually 
expressed as weighted integrals of the index of refraction structure constant, as described below. 
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Consider, for example, observation of a star with a large aperture telescope. As noted earlier, the large 
aperture negates the effect of scintillation by averaging over constructive and destructive interferences. 
Since light from a distant object can be well approximated by plane waves, the only prominent effect 
of the turbulent atmosphere is thus the change in the angle of arrival of the plane wave. The random 
arrival angle limits the resolution of the telescope. In fact, the resolution of the telescope is equivalent 
to the resolution of an aperture, ro (called the coherence length of the atmosphere), determined by the 
atmospheric structure constant and wavelength of light [5]: 


ro 


0.42fc 2 sec 0 


[ Z Cl(h)dh 

J ho 


- 3/5 


where k is the wave vector (2tt/X) of the plane wave; h is the height above sea level; and ho is the altitude 
of the location of the telescope. Prom basic principles of optics, an aperture with diameter ro has a 
full-width-half-maximum (FWHM) angular resolution of 


angular resolution = — 


For a wavelength of 0.5 jum, the coherence length ro is on the order of 10 cm. In experimental observations, 
it is convenient to define “seeing” as the FWHM of the angular spread of the star. For ro of 10 cm, the 
seeing is approximately 1 arcsec, or 5 /xrad. It must be noted that though subarcsecond seeing conditions 
are possible at TMF, a more representative value for the seeing at TMF is 2-3 arcsec, or an r 0 of a little 
under 5 cm. 

Given the AFGL CLEAR I model for the altitude dependence of the structure constant, the above 
equation for ro gives a single value for the atmospheric coherence length irrespective of the meteorological 
conditions. To better model the existing local weather conditions, we scale the structure constant function 
by a multiplicative constant so that the predicted value of ro matches the ro derived from measured 
seeing. The scaled C 2 is then used to determine other atmospheric-related quantities, such as scintillation 
variation and beam wander. For this purpose, the seeing was measured several times during a single run 
of the GOLD experiment. 

C. Beam Motion and Jitter 

To account for intensity fluctuations due to pointing inaccuracies and beam wander, we assume a 
Gaussian intensity profile for the uplink laser beam. That is, the normalized intensity pattern takes the 
form 


I{6 x ,6 y )=ex p 



/ fix + fiy \ 

V ©fed ) 


where S x and S y are the angular deviation or error in the pointing. If the transmitter and receiver are on- 
axis and there is no pointing offset (6 X = 6 y = 0), the receiver would see a power P 0 from the transmitter. 
Because of pointing jitter and atmospheric-induced beam motion, however, S x and fi y are random variables 
(RVs) and not constants. We assume fi x and fiy to be identically and normally distributed with zero mean 
and variance cr 2 . The variance cr 2 is the sum of the variance of atmospheric-induced and transmitter- 
induced pointing errors. The variance of the atmospheric-induced beam motion is again given by a 
weighted integral of the structure constant [5, Chapter 6], but as a first-order approximation, we assume 
it to be equal to the seeing. The transmitter-induced jitter is more difficult to ascertain, but is typically 
small relative to atmospheric effects. 
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The new random variable I can be shown to have the following cumulative probability distribution 
function (CDF) [6]: 


Pi (I < i) = F[(i) = i@ for 0 < i < 1 

where 



is a measure of the pointing accuracy with respect to the beamwidth. We immediately see that the mean 
and variance of I are given by 




P 

P+1 


and 


<I 2 > - <1 > 2 = 


P 

(P + 1) 2 


1 

P + 2 


respectively. It is obvious that a large value of p is desirable. In fact, for large values of p, the mean 
approaches the optimum value of 1 and the variance tends to 0 as 1/P 2 . For small values of p, (P « 1), 
on the other hand, both the mean and variance are proportional to p. 


To achieve a large value of P, either the angular beamwidth must be increased or the beam jitter 
must be decreased. With little control over the atmospheric-turbulence-induced beam wander, one is 
often limited in practice to increasing the beam divergence for obtaining a particular value of p. But 
the received power is inversely proportional to the square of the beam divergence. An optimum value 
of P, therefore, exists and depends on atmospheric conditions. The statistical analysis of the effect of a 
constant pointing offset is complicated by the loss of circular symmetry, and one may be forced to resort 
to numerical or Monte-Carlo simulations for a complete analysis. Thus, for simplicity, we will account 
for constant pointing offset only through an exponential loss factor, which is a reasonable assumption if 
the offset is much smaller than the beam divergence. 

D. Atmospheric Scintillation 

The pdf of the normalized received intensity, S , at a point receiver due to scintillation can take several 
forms. For weak turbulence, the pdf takes on a lognormal distribution. This can be understood from the 
fact that S is proportional to 


S oc |exp(x + i^)| 2 

where \ and <p are the amplitude and phase of the wave, respectively, and both are normally distributed 
random variables. The quantity 2y is often referred to as the log-irradiance. For strong turbulence over 
long paths, the intensity fluctuations are severe and the distribution is exponential. Andrews and Phillips 
have shown that the scintillation takes the form of an I-K (i.e., containing I and K Bessel functions) 
distribution, which reduces to the two limits given above for weak and strong turbulence [7,8], Since the 
GOLD experiments were performed at a relatively high altitude (2.2 km for TMF) and the source of most 
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turbulence is at or near sea level, we assume a weak turbulence model and use the lognormal distribution 
for S. That is, 


fs(s) 


1 

y/2ir of 


1 

- exp 
s 



(In s-l m ) 2 


where l m = — of/2 is the mean of the log-irradiance distribution. The S, being a normalized quantity, 
has a mean of 1. The variance of the log-irradiance is determined from the following weighted integral of 
the index of refraction structure constant [9,10]: 


erf = 2.24 k 7/6 (sec (9) 11/6 [ Cf(/i)/i 5/6 dh 

Jho 


where k, h, and ho are as defined before. Experiments have shown that erf does not arbitrarily increase 
with increasing atmospheric turbulence [11]. In fact, the scintillation variance reaches a saturation value 
of approximately 0.25of = 0.3 and even decreases for high levels of turbulence. The variance of S itself 
is given by 


= exp (of) - 1 


Note that it is customary to characterize the lognormal distribution of S using the mean and variance 
of the log-irradiance. Thus, neither the mean nor the variance of S is immediately apparent from /s(s). 
Given y = exp(x), where x is a Gaussian RV, and thus y is a lognormal RV, we use the following 
convention: By log-variance of y, we mean the variance of x. 

Combining the effect of atmospheric scintillation and beam-pointing errors, the received power can be 
written as 


Pr 


„ n^secO 

— V toeVtoeVo 


1 

* (ze bd ) 2 


We emphasize that in the above equation I and S are RVs. The pdf of the received signal resulting 
from the product of a lognormal and beta-distributed random variable can be derived analytically. This, 
however, involves the use of cumbersome integral functions. Through Monte-Carlo simulations, Kiasaleh 
has shown that, for a scintillation variance larger than 0.2, the pdf of the product of S and I can be 
approximated by a lognormal distribution [12]. We will, therefore, assume the received power Pr to be 
lognormally distributed to arrive at an analytical expression for the final result. 

The voltage readout, V, of the detection system on board the satellite is the product of the received 
power and the transimpedance gain, G, of the detector-amplifier combination. That is, V = G Pr. In 
principle, the statistics of the detector and amplifier must be taken into account for accurate analysis. 
For the case at hand, scintillation and jitter are dominant, and we therefore assume the statistics of V to 
be identical to Pr. The pdf of V, being lognormal, is 


fv{v) = 


1 

\/2'Ka 2 


1 

- exp 
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To determine the value of vq and a, we use the properties of independent RVs. That is, (V) = GPo(I)(S) 
and (V 2 ) = ( GPo) 2 {I 2 )(S 2 ). Moreover, from the properties of the lognormal distribution, In — —0.5 a 2 
+ ln(V) and a 2 = ln(V 2 ) — ln(V') 2 . From these expressions, we find that 


fj 2 = of + In 


(1 + / 3) 2 
( 2 + /?)/? 


and 


The photodetector/amplifier gain factor G has units of volts per watt. The mean and variance of the 
received signal are not independent, as they both depend on P and of. For example, an increase in 
atmospheric turbulence resulting in a larger scintillation variance will decrease the mean. 

E. Effect of Multiple Uplink Beams 

The adverse effect of scintillation can be reduced by the use of multiple uplink beams, each incoherent 
with the rest and separated by a distance larger than the atmospheric coherence diameter, rQ. For such a 
case, the signal arriving at the detector is the sum of the signal from each uplink beam with independent 
identically distributed (i.i.d.) statistics. That is, 


1 N 

s = jvE s - 


71=1 


where N is the number of independent beams, and the random variables S n are lognormally distributed. 
We assume that the total uplink power remains constant and, thus, each beam has 1/Nth of total power. 
From the properties of sums of i.i.d. random variables, we find that S' is a convolution of N lognormal 
distributions. Analytic expressions for the convolution of lognormal functions are unavailable, and thus we 
resort to numerical techniques using characteristic functions and fast Fourier transforms (FFTs). Note 
that for a large N, according to the central limit theorem (CLT), the distribution of S approaches a 
Gaussian. 

To illustrate the benefits of launching multiple beams during uplink, we choose, as an example, a mean 
of 100/ AT and a log- variance of 0.6 for each S n . Figure 3 shows the expected pdf when total laser power 
is equally distributed in 1, 2, 4, 8, or 16 beams. Several features are worth noting. First, the means for 
all the distributions are the same. Second, the location of the peak of the pdf approaches the mean of S 
when N is increased. Finally, the pdf is concentrated around the mean for large values of N. In short, 
though the mean varies little with an increasing number of beams, the variance drops significantly with 
additional beams. In fact, for the example in the figure, the standard deviations are 88, 64, 45, 32, and 
23 for 1, 2, 4, 8, and 16 beams, respectively. That is, the standard deviation is approximately one-half 
with 4 beams and one-quarter with 16 beams. 

In the presence of strong scintillation, dividing the laser power into several beams is essential for 
avoiding deep fades and surges. The number of beams needed to achieve a given bit-error rate (BER) will 
depend on the strength of the scintillation. One can, in principle, increase the number of beams arbitrarily 
to improve signal statistics at the receiver. Practical considerations, such as efficiency of dividing the laser 
power, modulation timing accuracy, and, most importantly, the available space, however, will typically 
limit the number of beams to less than 10. 
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Fig. 3. The pdf of the received signal intensity, S, showing the 
improvement in variance with multiple uplink beams. 


Part of the objective of GOLD is to understand the effects of multibeam propagation through the 
atmosphere. In the GOLD experiments, light from the laser is split into two and the path length of one 
of the beams is increased by a distance greater than the coherence length of the laser. The two beams are 
then launched into the coude so as to emerge from the primary separated by a distance of 20 cm (which 
is much greater than r 0 ). Thus, the actual variance of the received signal will be less than that predicted 
by the single-beam model. 


III. Link Budget 

Table 1 shows the link budget for transmission of a single high-power CW beam from TMF up to the 
satellite on November 17, 1995. The link budget is for the quad-photodetector (QPD) on the satellite. 
All data concerning the laser communication package on the satellite are provided by the Communication 
Research Laboratory (CRL) in Japan. Parameters pertaining to the TMF transmitter were determined 
experimentally. Atmospheric-related parameters were estimated from “seeing” conditions, as described 
earlier. We assume that the beam wander is dominated by atmospheric-induced motion and that it is 
equal to the seeing. Finally, background radiation from Earth reaching the detector is negligible, as the 
experiments are conducted at night from TMF. Values of parameters such as range, elevation, and seeing 
represent average values for the duration of the experiment. The atmospheric scintillation contribution 
to the log-variance is 0.52 in the above calculation. The remainder of the contribution to the variance is 
from atmosphere-induced beam wander. The effect of two beams on the expected variance is discussed 
in the next section. 


IV. Experimental Data and Analysis 

Figure 4 shows the QPD signal versus time for an approximately 100-min time period when two CW 
beams (each with ~6.6 W of power) were sent up to the satellite. The data were collected on the 17th of 
November 1995. The QPD signal was sampled every second and sent down through RF telemetry. Fur- 
thermore, the QPD has a gain value of 0.22 V/nW in the “low” gain setting, a gain value of 4.60 V/nW 
in the “high” gain setting, and few other intermediate settings. During the time corresponding to Fig. 4, 
the gain was set to low. Figure 5 shows the histogram of the signal shown in Fig. 4. The most likely 
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Table 1. Link budget for a single beacon uplink from TMF to ETS-VI. 


Parameter 

Value 

Equivalent value, 
dB/dBm 

Range (satellite to TMF distance) 

37,850 km 

— 

Zenith angle 

43.8 deg 

— 

Seeing 

6 prad 

— 

Transmitter information 



Laser wavelength 

514.5 nm 

— 

Laser power 

13,200 mW 

41.21 dBm 

Optical transmission 3 

0.75 

-2.22 dB 

Beam divergence 

20 /irad 

— 

Pointing loss 

4 fj , rad 

-0.09 dB 

Pointing jitter 

0 /rrad 

— 

Free space/atmosphere 



Propagation loss 

— 

-56.53 dB 

Estimated atmospheric transmission at zenith 

0.80 

— 

Atmospheric loss (at zenith angle) 

— 

-1.34 dB 

Receiver information 



Aperture diameter 

7.50-cm 

-23.55 dB 

Optical transmission 

0.15 

-8.24 dB 

Received power (no turbulence) 

10.49 nW 

-49.79 dBm 

Sensitivity (“high” gain) 

631.0 pW 

—62.00 dBm 

Margin 

16.62 

12.21 dB 

Reduction due to scintillation 

0.77 

-1.11 dB 

Reduction due to jitter 

0.71 

-1.48 dB 

New received power 

5.7 nW 

—52.38 dBm 

New margin 

9.0 

8.89 dB 

Net log-variance, a 2 

0.69 

— 

Probability of detection with above sensitivity 

99.6 percent 

— 

Scintillation fades greater than 3 dB 

33.7 percent 

— 

a Estimate based on losses in optical train components and reflectivity of 
telescope mirrors. 


value of the QPD signal was approximately 1.2 V. The clustering of points around the maximum value 
of 4.5 V indicates saturation of the detector due to scintillation-induced surges, while clustering near 
the minimum value of 0 V indicates the beam was outside the field of view of the QPD due to tracking 
errors. The fact that the statistics of the received signal are not lognormal but rather a convolution of two 
lognormal distributions makes analyzing the data difficult because of the lack of an analytic expression 
for the convolution. The histogram data may nonetheless be fit to a lognormal function to obtain an 
approximate value for the mean and log-variance. To do so, the histogram data were fit to a scaled 
lognormal distribution, A f(v), in an indirect way. The lognormal distribution can be written as a 
second-order polynomial by appropriate transformation of variables. Let w = In f(v) and u = \nv. By 
taking the log of the lognormal distribution function, we have 

w = au 2 +bu + c 
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Fig. 4. The QPD signal reading versus time in minutes since the start of the 
experiment (10:40 UTC) on November 17, 1995. 



Fig. 5. Histograms of the QPD signal reading. The asterisks show 
the number of times a given signal level was observed in the time 
interval of interest. The solid line is a fit of the histogram data to a 
log-normal function. 

The amplitude, mean, and variance of the lognormal distribution can be inferred from the parameters 
a, b , and c in a straightforward manner, the results of which are given below: 

1 

<7 = T== 
v— -2a 

lnz/Q = cr 2 (l -I- b) 


, , (Ini'o) 2 1 , 

In A = c+ ^-^2 h — ln(2-7r<7 2 ) 

Using least-square methods, the log of the raw data can, therefore, be easily fit to a polynomial to 
determine the mean and variance of the lognormal distribution. This technique eliminates the need 
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to take several partial derivatives of the lognormal distribution and simplifies the fitting algorithm. 
The fit is shown along with the original data in Fig. 5. Data points near saturation and data points 
around the minimum signal level resulting from pointing errors were not considered. The resul ting 
fit had a mean of uq = 1.81 V and a 2 = 0.49. These numbers account for contributions from both scintil- 
lation and beam motion/jitter. In the presence of multiple beams, it is more convenient to use the actual 
mean and standard deviation of the distribution as opposed to indirect measures such as log-variance. 
From In < V >= 0.5cr 2 + lnuo, we find that the mean of the fit is < V >= 2.3 V. It is also useful to 
define a normalized standard deviation, SD = (< V 2 > - < V > 2 ) 1 / 2 / < V > for the distribution. It is 
easy to show that SD = (exp(er 2 ) — l) 1 / 2 , and therefore the standard deviation is 0.80 for the fit. 


V. Comparison of Theory and Experiment 

Figure 6 shows the predicted (with one and two beams) and experimentally observed pdf’s for the QPD 
signal. Figure 6 also lists the mean and standard deviation for each of the theoretical and experimental 
curves. For two beams, the theory predicts the standard deviation of the signal to be 0.71, while it is 0.80 
for the experimental fit. The difference arises mainly from the uncertainty in the transmitter pointing 
jitter and the unaccounted satellite tracking errors (discussed later). The electronics on board the satellite 
are susceptible to Earth’s radiation belts, and the resulting noise can increase signal variation as well. 
Additional data from the planned experiments with a single uplink beam will enable better comparison 
to theory. 



Fig. 6. Expected and observed pdf of the QPD signal, showing the 
experimental lognormal fit of Fig. 5, the theoretical pdf for a single 
uplink beam, and the expected pdf in the presence of two uplink 
beams. 


The discrepancy between the predicted and observed mean values, on the other hand, is just over 
1.5 dB. It is interesting to note that the theoretical and observed most-likely value of the QPD voltage 
reading is almost identical. The agreement between theory and experiment for the mean value is surpris- 
ingly good given the various unknown factors that are currently being investigated. Recent experiments 
suggest a reduction in the sensitivity of the detector with time due to temperature increases. The effect of 
excessive exposure to Earth’s radiation belts on the QPD and charge-coupled device (CCD) sensitivity is 




still unclear. Some of the difference can be attributed to uncertainties in the optical loss of the transmitter 
system, which was estimated from the number of elements in the optical train. Also, the atmospheric 
transmission value at the uplink wavelength needs to be updated from the planned recalibration of the 
atmospheric visibility monitoring (AVM) data. Moreover, neither the transmitter pointing offset nor the 
beam divergence is accurately known. 

As is apparent from Fig. 4, there is a slow (on the order of minutes) decay in the signal level followed 
by a sharp jump to a high value. This is due to drift in the tracking of the beacon by the satellite, which 
is not included in the analysis. The tracking errors cause the out-of-focus beam to fall partially outside 
the area of the QPD and, hence, cause a reduction in the detected signal. The slow decays in received 
signal power partly explain the large variance in the observed signal. Faster sampling of the beacon signal 
from the communication detector should provide additional data about temporal statistics of scintillation. 
Downlink signal statistics will provide information on satellite vibrations, pointing accuracy, and pointing 
drift that will further enable us to better isolate atmospheric scintillation effects. 
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An array feed combining system for the recovery of signal-to-noise ratio (SNR) 
loss due to antenna reflector deformation has been implemented and is currently be- 
ing evaluated on the Jet Propulsion Laboratory 34-m DSS-13 antenna. The current 
signal-combining system operates under the assumption that the white Gaussian 
noise processes in the received signals from different array elements are mutually 
uncorrelated. However, experimental data at DSS 13 indicate that these noise pro- 
cesses are indeed mutually correlated. The objective of this work is to develop a 
signal-combining system optimized to account for the mutual correlations between 
these noise processes. The set of optimum combining weight coefficients that maxi- 
mizes the combined signal SNR in the correlated noises environment is determined. 
These optimum weights depend on unknown signal and noise covariance parameters. 
A maximum-likelihood approach is developed to estimate these unknown parame- 
ters to obtain estimates of the optimum weight coefficients based on residual carrier 
signal samples. The actual combined signal SNR using the estimated weight coef- 
ficients is derived and shown to converge to the maximum achievable SNR as the 
number of signal samples increases. These results are also verified by simulation. 
A numerical example shows a significant improvement in SNR performance can be 
obtained, especially when the amount of correlation increases. 


I. Introduction 

An array feed-combining system has been proposed for the recovery of signal-to-noise ratio (SNR) 
losses caused by large antenna reflector deformations at Ka-band (32-GHz) frequencies in the Deep 
Space Network [1]. In this system, a focal plane feed array is used to collect the defocused signal fields 
that result from these deformations. All the signal power captured by the feed array is then recovered 
using real-time signal-processing and signal-combining techniques. The optimum combiner weights that 
maximize the combined signal SNR were derived in [1] under the assumption that the white Gaussian 
noise processes in the received signals from different array elements are mutually uncorrelated. These 
optimum weights depend on unknown signal and noise parameters that need to be estimated. The work in 
[1] proposed to estimate the optimum weights from the observed residual carrier received-signal samples 
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using a maximum-likelihood (ML) approach. The actual combined-signal SNR in an uncorrelated noises 
environment when the estimated weights are used in place of the optimum weight coefficients was also 
derived in [1]. 

A seven-element array feed combiner system is currently being evaluated at the JPL DSS-13 34-m 
antenna. Although the work in [1] assumed mutually uncorrelated noise processes, experimental data [2] 
indicate that the noise processes in the received signals from different feed elements are indeed correlated, 
with correlation coefficients of the order of 0.01 under clear sky conditions. Since the noise in each of 
the array feed element signals consists of receiver white noise plus noise due to background radiation, it 
has been conjectured that this small correlation is caused by near-field atmospheric background noise. 
Recent data gathered at DSS 13 in more adverse weather conditions, however, indicate both increases 
and decreases in the observed amount of correlation. This experimental work is still in progress. Larger 
correlations may also result from undesired radiation source emissions gathered by the antenna side lobes. 
Moreover, our recent work [5] has shown that the array feed combining system derived in [1], which is 
suboptimal in the correlated noises environment, actually can have a better performance in the presence 
of correlated noises. Therefore, it is important to develop optimum signal-combining techniques that 
account for the mutual correlation between the noises in the signals from different array elements. That 
is the objective of this work. 

As a first step towards this objective, in Section II we provide a derivation of the set of optimum 
combining weight coefficients that maximizes the combined-signal SNR in this correlated noises environ- 
ment. These optimum weights depend on unknown signal amplitude and phase parameters as well as 
noise variance and correlation parameters. An ML approach is then developed in Section III to estimate 
these unknown parameters and arrive at an ML estimate of the optimum weight coefficients based on 
residual carrier received-signal samples. The actual combined-signal SNR is derived in Section IV when 
the estimated weights are used in place of the optimum weights, with the details given in Appendices A 
and B. The SNR performance is shown to converge to the maximum achievable SNR as the number of 
signal samples used in the estimates increases. These results are also verified by simulation. Numerical 
examples are given in Section V with a particular choice of noise covariance matrix and signal parame- 
ters and show a significant improvement in SNR performance compared to the previous signal-combining 
system developed in [1], especially when the amount of correlation increases. 


II. Array Feed Signals and Optimal Combining Weights 

Consider a A-element array and the NASA Deep Space Network standard residual carrier modulation 
with binary phase shift keyed (PSK) modulated square- wave subcarrier [4]. The received signal from each 
array element is downconverted to baseband and sampled. Similar to the combining system proposed in 
[1], only the residual carrier portion of the received signal spectrum will be used to estimate the unknown 
parameters in the combiner weights. The full-spectrum modulated signals from the array elements, which 
contain both the modulated sidebands as well as the residual carrier spectrum, are subsequently combined. 
As in [1], assume that the higher bandwidth primitive baseband signal samples are lowpass filtered by 
averaging successive blocks of M B samples to yield a full-spectrum signal stream B for each array element. 
Additive white Gaussian noise is assumed to be present in the primitive baseband signal sequences from 
each of the array elements. The white Gaussian noises in the primitive baseband samples from different 
array elements are assumed to be mutually correlated. Specifically, the noise samples corresponding to 
different array elements are assumed to be mutually correlated at any given time instant, but uncorrelated 
at different time instants. Let 

Vk{iB) = V k [cos 6 + js(i B ) sin 6} +n k (i B ), i B = 1,2,- (1) 

denote the stream B signal samples from the fcth array element. The complex signal parameters Vfe, 1 < 
k < K, represent the unknown signal amplitude and phase parameters induced by the antenna reflector 
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deformation. Moreover, 5 is the modulation index, s(i B ) = ±1 is the transmitted data, and {n k (i B )} 
is the zero-mean white Gaussian noise corruption in the stream B signal samples from the fcth array 
element. The primitive baseband signal samples are also more narrowly lowpass filtered by averaging 
successive blocks of Ma samples to yield a residual carrier signal stream A for each array element. 
Clearly, Ma > Mb and r? = Ma/Mb is the ratio of the bandwidth of stream B to stream A. Let 

u k (i A ) =V k cos 6 + m k (i A ), ^a — 1,2,--- (2) 

denote the stream A signal samples from the A;th array element. Here is the zero-mean white 

Gaussian noise corruption in the stream A signal samples from the fcth array element. Let A*, A 1 ' , 
and A t denote the complex conjugate, the transpose, and the complex conjugate transpose of the 
matrix A, respectively. In order to specify the correlations between the white noise sequences cor- 
responding to different array elements, consider the noise vectors n(i B ) = (ni(i B ) • • • nx(i B )) T and 
UBJAa) — (^i(*a) • • • Then {n(i B )} and {to(m)} are each sequences of independent identi- 

cally distributed (i.i.d.) zero-mean complex Gaussian random vectors of dimension K. The respective 
covariance matrices R B = {r B kj} = E[n(is)n(is) t ] and R A = {rAkj} = E[m(i / i)m(* J 4 ) t ] of n(i B ) and 
m(iA ) specify the mutual correlations between the white noises in the signal streams from different array 
elements. For example, r Bk j is the correlation between the noise variables n k (i B ) and rij(i B ). Because 
of the different averaging rates in streams A and B on the primitive baseband signals, it follows that 
R b = t\R a - Finally, these different averaging rates also imply that m(M) is independent of n(i B ) pro- 
vided that ia < i-B and the samples averaged to yield m(iA ) occurred prior to the samples averaged to 
yield n(i B ). 

Application of the complex combining- weight coefficients i W k , 1 < k < K yields the combiner output 
sequence z c (i B ) — s c (i B ) + n c (i B ), where 


K 

Sc (i B ) = J2 w ^ eMiB)s 

k = 1 


(3) 


and 


K 

n c (i B ) = 'Y^W k n k {i B ) 
k = 1 


(4) 


are the signal and noise components, respectively. Define W = ( W\ ■ ■ ■ Wk) f to be the vector of 
combining-weight coefficients. The objective is to determine the optimum weight vector W that maxi- 
mizes the SNR of the combiner output defined by 


, W \ = l E Mfg )]| 2 lsc(fi ?)| 2 
Var[z c (i B )] Var[n c (i B )] 


(5) 


The optimum weight vector and the maximum achievable SNR have been derived previously in [3]. For 
the sake of completeness, we provide a derivation below that is slightly different from that in [3] . Define 
V — (Vj • • • Vk) T to be the vector of complex signal parameters. Then, from Eqs. (3) and (4), we have 
|s c (*b )| 2 = \W T V\ 2 and Var[n c (i B )] = W t R f W*. substitution of which in Eq. (5) yields 


7(00 = 


\W T V\ 2 

W t R b W* 


( 6 ) 
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Since R B is a positive definite Hermitian matrix, there is a unitary matrix Q such that 

Rb = $D 2 Q (7) 

where D is a real-valued K x K diagonal matrix with the fcth diagonal term given by y/r B kk- Using 
Eq. (7), we have 


W t R b W* = (w T Q^Dj ( DQW *) 

= (DQW*y (DQW*) = \\DQW *\\ 2 


Moreover, since Q is the inverse of the matrix , 


(m*) 


T' 


-1 



= D~ l Q 


Hence, by using Eq. (9), we can write 


W T V = W T (DQ*f \(DQ*y 




= (DQW*y (D~ l QV) 


( 8 ) 


(9) 


( 10 ) 


So substituting Eqs. (8) and (10) in Eq. (6) and applying the Schwartz inequality gets 


7(10 = 


I (DQW*y (D-'QV) I 2 
\\DQW *\\ 2 


WDQW-Tm^QVW 2 

\\DQW *\\ 2 


= \\Dr l QV \\ 2 =1MAX 


( 11 ) 


Moreover, equality holds in Eq. (11) if and only if for some complex-valued constant a, DQW* = 
a (D~ 1 QVy . So, by using Eq. (7), the set of optimum weight vectors W OPT that achieves 7 (W OPT ) = 
Imax is given by 


Wqpt = a ( &DT 2 Q ) V - a (R* B )~ l V* = ^ CSk) -1 1* (12) 

where a is an arbitrary complex-valued constant. Moreover, it follows from Eqs. (7) and (11) that the 
optimum SNR 7 max can be written as 
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1MAX = (D-'QVy (D-'QV) = V?R b 1 V 


(13) 


Note that, in the uncorrelated noises case, R B l — D 2 . So the set of optimum weight vectors, Eq. (12), 
in this case is given by W OPT — aD_~ 2 V*, and the optimum SNR is given by 

7 MAX = WD-'VW 2 = E {r^ = 7 (14) 

tZ i rBkk 

which is the sum of the array element output SNRs. These results for the uncorrelated noises case agree 
with previous results derived in [1]. 


III. Parameter Estimation 

The signal parameter vector V_ and the noise covariance matrix R A are not known and need to be 
estimated to obtain an estimate of one of the optimum weight vectors, W OPT . given by Eq. (12). Assuming 
that these unknown parameters are not random, we propose to use ML estimates based on the stream A 
residual carrier-signal vector samples {u(iji)}, where u(z^) = (u\(i A ) ■ ■ ■ UK(iA)) T ■ Instead of estimating 
V directly, consider estimating X = VcosS. Note from Eq. (2) that (u('m)} is an i.i.d. sequence 
of complex Gaussian random vectors with mean X and covariance matrix R A . It then follows from 
multivariate statistical analysis [6,7] that, based on observations {u(iA — 1), • • • ,u(i A — L)}, 

x ML {i A ) = \ £ a(0 ( 15 ) 

l=iA ~ L 

and 

^A,Mi(^) = ^ E [n(0 — ^Ml(^)] [n(0 _ =^Ml(^)] (16) 

1—lA — -k 

are the respective ML estimates of X and R A . By the invariant property of ML estimators, 

Y-ml^a) = ^-^Kml^a) (17) 

is then the ML estimate of V and 

W-opt^a) = — (E-aml^a)} Vml^a) (18) 

is the ML estimate of W OPT given by Eq. (12). 

It is well known [6,7] that X_ Mpj (i A ) is a complex Gaussian random vector with mean X_ and covariance 
matrix (1 /L)R A . Furthermore, X_ml(^a) and R mb (ia) are statistically independent, and LR ML {i J \) has 
the same distribution as the random matrix 
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( 19 ) 


i= 1 

where Z_ i is a sequence of i.i.d. zero-mean complex Gaussian random vectors with covariance matrix R A . 
This type of distribution is called a complex Wishart distribution with parameters R A and (L — 1), and A 
in Eq. (19) is said to have a CW(R A , L — 1) distribution [7]. It has been shown in [8] that, for L > K + 1, 

E W-‘I - « 


for a j K x K CW(R A ,L — 1) distributed random matrix A. Since (A * ) ~ 1 = (A -1 )*, it then follows from 
using the property of Eq. (20) for the complex Wishart matrix LR ML {i A ) along with Eqs. (12), (17), and 
(18) that, for L > K + 1, 


E 








olL 


V(L-K-I) 


R A ~ l V* 


L 

L-K-l 


W-OPT 


(21) 


Hence, the ML estimate W OPT (i , if of the optimum weight vector W OPT is actually a biased estimate. 
This is not a problem, since it is clear from Eq. (12) that the optimum weight vectors are not unique and 
any complex scaled version of an optimum weight vector is also optimum. So the constant a in Eq. (18) 
can be set arbitrarily. We shall set a = (L — K — l)/L for the purpose of normalization and also assume 
that L > K + 1. The ML estimate of the optimum weight vector that will be used here is, therefore, 


with mean 


Kml^a) 


L-K-l 

Lr) 

L-K-l 
Lr] cos 6 


(&A ,ml (U))~ Vul^a) 
(J^a,ml (M))" 





( 22 ) 


(23) 


that is also an optimum weight vector. 

Note that V_ ml (ia) and Ra,m l(m) are both consistent estimates, i.e., they both converge with prob- 
ability one to their respective expected values V and R A in the limit as L tends to infinity. So it follows 
that W_ ML (i A ) also converges with probability one to the optimum weight vector W Q in the limit as the 
number L of samples tends to infinity. The analysis in Appendix C shows that this convergence also holds 
in the mean-square sense. These properties indicate that we may expect the actual combiner output SNR 
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using the estimated weight vector W ML (i a) to " also converge to the maximum achievable SNR 7 max as 
L tends to infinity. That result will be shown in the next section, which derives an explicit expression for 
the actual combiner output SNR for finite L. 

As in [ 1 ], these weight coefficient estimates are used in a sliding- window structure to produce the 
following combiner output sequence: 


„ T 

z c {iB) = W ML (i A )y(i B ) 


(24) 


where i A is the largest integer less than i B , so that the residual carrier-signal vector samples 
{u (i A — l) , • • ■ , u (i A — L) } used for estimating W ML ( i A ) occur before the full-spectrum signal vector 
sample y(i B ) = {yi{i B ), • ■ • , y/c(*e)) • This ensures that the noise vector n(i B ) in y(i B ) is statistically 
independent of {u (i A — l) , ■ • • ,u(i A — L) }, and hence the statistical independence between W ML (i^) 
and n(i B ). 


IV. SNR Performance Analysis 

The combiner output, Eq. (24), can be written as 


z c (* b ) = s c (i B ) + n c (i B ) 


(25) 


where 


s c (i B ) = W T ML (i A )Ve^ s 


(26) 


and 


n c {i B ) = W_ ML (i A ) n(i B ) 


(27) 


are the signal and noise components, respectively. In the following analysis, V , R B , and s(i B ) are assumed 
to be nonrandom parameters. As noted above, W ML (i A ) and n(i B ) are statistically independent. Since 
the components of n(i B ) are all of zero mean, it follows from Eqs. (26) and (27) that n c (i B ) also has zero 
mean and moreover is uncorrelated with s c (i B ). Thus, it follows from Eq. (25) that the actual SNR 7 ml 
of the combiner output can be written as 


7 ML 


1 E[*c (is)] I 
Var[z(i B )] 


1 e M^b)] I 

Var[s c (j B )] + Var[n c (iB)] 


Now, it follows from Eqs. (26) and (23) that 


E[s c M] | 2 = | E [W T ML (i A )]V | 2 = | WlV | 2 = | V}R~ b x V I 2 = r MAX 


(28) 


(29) 


where 7 max is the maximum achievable SNR given by Eq. (13). The explicit expression, Eq. (A-9), 
for Var[s c (i B )] is derived in Appendix A. Moreover, the derivation in Appendix B yields the expression, 
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Eq. (B-5), for Var [n c (i B )}. So, by using Eqs. (29), (A-9), and (B-5) in Eq. (28), one arrives at the 
following expression for the actual combiner output SNR performance: 


7 ML 


7 MAX 

1MAX{1 + Cl) + C-i 


(30) 


where 


LK — K 2 — K + 1 (L- A- !)(£-!) f 1 \ 1 

1 (L - K)(L - K - 2) (L - K)(L -K- 2) V^cos 2 6 ) L-K- 2 1MAX 

K(L-l)(L-K-l) 1 

2 L(L - K){L - K - 2) r]cos 2 6 K 1 

Since both C\ and Ci converge to zero as L tends to infinity, it follows from Eq. (30) that the actual 
combiner output SNR 7 ml converges to the maximum possible achievable SNR, 7 max, as the number 
of signal samples, L, used in the estimates tends to infinity. 


V. Numerical Example 

We consider here a numerical example using a K = 7 element array feed. In this example, a modulation 
index, S = 80 deg, and a primitive sample period of To = 2.5 x 10~ 8 s are assumed. The full-spectrum 
modulation signal is assumed to be of bandwidth 2 x 10 6 Hz, which yields M B = 20. Moreover, the 
ratio of the full-spectrum bandwidth to the residual-carrier bandwidth is 77 = Ma/Mb = 200. A nominal 
P T /N 0 of 65 dB-Hz is considered with a corresponding 7 = ( Pt/N 0 )MbT 0 (recall that 7, which is given by 
Eq. (14), is the sum of the array element SNRs). The white Gaussian noise processes in the received signals 
from different array feed elements are assumed to be correlated. Moreover, the correlation magnitudes 
are assumed to be inversely proportional to the distances between feed centers (this assumption has not 
been verified for the array feed system at DSS 13). In the numerical examples below, the following noise 
covariance matrix is considered: 

H B = 


- 

1 
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(32) 


Figure 1 shows the K = 7 element feed array geometry and the relative distances between feed centers. 
The main feed is labeled feed 1 and is surrounded by the six others. In the correlation matrix, Eq. (32), 
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Fig. 1. Feed array geometry 
for K = 7 feeds. 


the noise power in the feed-element stream B received signal samples are all assumed to be equal and 
normalized to one ( r B kk — 1 f° r all k). The maximum possible noise correlation (and also correlation 
coefficient) magnitude is denoted by p max ■ The corresponding noise correlation magnitudes in this matrix 
are identical and reflect an inverse dependence on distances between feed centers. The noise correlation 
magnitudes between nearest neighbor feed pairs are equal to the maximum correlation magnitude p max 
(for example, between the main feed, feed 1, and any of the outer feeds). As can be seen from Fig. 1, 
the next nearest neighbor feed pairs (for example, feed 2 and feed 4) are of a distance equal to \/3 times 
the distance between nearest neighbor feed pairs, and hence have a noise correlation magnitude equal to 
Pmax/VS- Finally, the furthest feed pairs (for example, feed 2 and feed 5) are of a distance equal to twice 
the distance between closest feed pairs, and so have a noise correlation magnitude equal to p max /2. The 
parameter <j> in Eq. (32) specifies the noise correlation phases. The rationale for assigning these phases 
will be described below. Since the noise power in each feed element is equal to one, the sum of the feed 
element SNRs is equal to the total received power, and so 7 = Pt- For the complex signal parameters 
vector, F, we shall use 


y : 


a/ 


' 6 


■j&l L 

e J * 





6 


e 3 



(33) 


where 0 < (3 < 1 represents the fraction of total received signal power in the main antenna feed (feed 1) 
as defined in [9]. The remaining total received signal power is then evenly distributed among the other 
six feed elements (this assumption does not appear to be generally valid for severely distorted antenna 
reflectors). The signal phases in Eq. (33) were chosen arbitrarily. The parameter <p in Eq. (32) specifies 
the relation between the noise correlation phases and the signal parameter phases as follows: Let Vk 
denote the fcth component of the vector V given by Eq. (33). Then for k > j, <j> = [phase of rskj] 

— [(phase of Vk) — (phase of V))]. Moreover, since R B is a Hemitian matrix, —<j> = [phase of r B kj] 

— [(phase of Vk) — (phase of V})] when j > k. That is, (j) represents the difference between the phase of 
the noise correlations and the phase differences between the corresponding signal components in V. The 
significance of this will be discussed below. 

In the numerical examples below, the signal combining scheme given by Eq. (24) will be referred 
to as the correlated noises algorithm. The signal-combining scheme from [1] will, on the other hand, be 
referred to as the uncorrelated noises algorithm. A comparison of the SNR performance of these two signal 
combining schemes will be made below, assuming the noise covariance matrix given by Eq. (32) and the 
signal parameter vector given by Eq. (33). Our previous work in [5] determined the SNR performance of 
the uncorrelated noises algorithm in the environment where the noise processes in the feed element received 
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signals are indeed correlated. For a given signal vector V and set of noise correlation magnitudes, the 
results in [5] show that the best-case performance of the uncorrelated noises algorithm occurs when each 
of the noise correlations between feed pairs has a phase that is exactly 180 deg from the phase difference 
between the corresponding signal components in V. For the correlation matrix, Eq. (32), this “best case” 
scenario for the uncorrelated noises algorithm occurs when the parameter (p — 180 deg. Moreover, the 
worst-case performance was shown in [5] to occur when each of the noise correlations between feed pairs 
has a phase that is exactly equal to the phase difference between the corresponding signal components in 
V. This “worst-case” scenario for the uncorrelated noises algorithm occurs when the parameter cj> = 0 deg 
for the correlation matrix, Eq. (32). The numerical examples below will consider these two extreme cases 
only. The performance comparisons between these two algorithms can then be made under both the 
most-favorable and the least-favorable situations for the uncorrelated noises algorithm. 

The SNR performances will be compared in terms of the combining gain, which is the ratio of the 
actual SNR performance to the sum of the SNRs of the array feed element received signals ( 7 ml /7 
for the correlated noises algorithm). Note that 7 is also the maximum achievable SNR in the uncor- 
related noises environment. Hence, the combining gain represents the SNR gain relative to the best 
possible performance in the uncorrelated noises environment. We shall, therefore, refer to an SNR 
gain when the combining gain is positive (in dB) and to an SNR loss otherwise. The SNR perfor- 
mance analysis given in Section IV showed that the combining gain 7ml/7 of the correlated noises 
algorithm converges to the maximum possible combining gain, 7 max /l, in the limit as the number of 
samples, L , approaches infinity. Figure 2 plots this maximum possible combining gain as the phase 
parameter (p in Eq. (32) is varied between 0 and 180 deg for values of p max equal to 0.1, 0.15, and 
0.2, and when (3 = 0.7. These results show that the best- and worst-case scenarios for the correlated 
noises algorithm are the same as those for the uncorrelated noises algorithm. These results show that 



Fig. 2. Maximum possible combining gain versus with /? = 0.7 
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the maximum possible combining gain is positive for a majority of the <j> phase values. Moreover, the 
maximum possible combining gains increase with increasing p max . Simulations were also performed to 
validate the analytical result, Eq. (30), that yields the combining gain for the correlated noises algorithm. 
Table 1 compares the simulated combining gain to the analytical result for the best-case R B matrix (with 
<t> — 180 deg) at various values of p max and L. Table 2 displays the corresponding comparisons for the 
worst-case R B (with cf> — 0 deg). The simulated combining gain from these two tables can be seen to be 
within 3 percent of the analytical results. This appears to validate the analysis. 


Table 1. Correlated algorithm simulated 
combining gain for best-case R B . 


SNR Gain, dB 

L 

Pmax 

Simulation 

Analytical 

500 

0.050 

0.85896 

0.80340 

500 

0.100 

2.24048 

2.17891 

500 

0.150 

4.75655 

4.67706 

500 

0.175 

7.25759 

7.14633 

500 

0.200 

15.06133 

14.60054 

5000 

0.050 

0.88286 

0.87982 

5000 

0.100 

2.26456 

2.25971 

5000 

0.150 

4.78108 

4.77260 

5000 

0.175 

7.28267 

7.26961 

5000 

0.200 

15.08713 

15.03570 

10 6 

0.050 

0.86770 

0.88821 

10 6 

0.100 

2.24674 

2.26859 

10 6 

0.150 

4.76168 

4.78312 

10 6 

0.175 

7.26309 

7.28324 

10 6 

0.200 

15.06514 

15.08578 


In Fig. 3, the combining gains of both algorithms are plotted versus the number of samples, L , for 
(3 = 0.7 (70 percent of the power in the main feed), using the best-case R B and with p mox having values of 
0.05, 0.1, and 0.15. We notice that for a fixed pmax » the correlated noises algorithm (CNA in the figures) 
actually has a smaller combining gain than the uncorrelated noises algorithm (UNA in the figures) for 
small values of L below a threshold value. In the examples considered, this threshold value increases 
with decreasing p max and is always less than about L — 500. Note that although the ML estimator is 
asymptotically optimal as L — » oo, it is not necessarily the best estimator for small values of L. Hence, 
the correlated noises algorithm may not have the best possible performance at small values of L. The 
convergence of the combining gain to the optimum theoretically achievable value was proved in the last 
section. These examples show that convergence to within 0.1 dB of the limiting value occurs at about 
L = 1000 samples. Figure 4 shows the respective combining gains as a function of f3 , the fraction of total 
received signal power in the main feed, for L = 5000 samples and values of pmax equal to 0.05, 0.1, 0.15, 
and 0.2. We can see that the performance of the correlated noises algorithm is much less sensitive to (3 in 
this example than that of the uncorrelated noises algorithm. The performance superiority of the corre- 
lated noises algorithm also increases significantly with increasing (3. The respective combining gains versus 
Pmax with L = 5000 and (3 taking on values of 0.7, 0.8, and 0.9 are shown in Fig. 5. This figure shows that 
the performance improvement of the correlated noises algorithm over the uncorrelated noises algorithm 
increases significantly with increasing p max - The results in these three figures indicate that there can be 
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Table 2. Correlated algorithm simulated 
combining gain for worst-case fig. 


SNR Gain, dB 

L 

Pmax 

Simulation 

Analytical 

500 

0.050 

-0.63659 

-0.68802 

500 

0.100 

-1.05467 

-1.10542 

500 

0.150 

-1.33373 

-1.38419 

500 

0.175 

-1.43041 

-1.48083 

500 

0.200 

-1.50179 

-1.55220 

5000 

0.050 

-0.61306 

-0.61410 

5000 

0.100 

-1.03140 

-1.03182 

5000 

0.150 

-1.31079 

-1.31072 

5000 

0.175 

-1.40765 

-1.40738 

5000 

0.200 

-1.47920 

-1.47877 

10 6 

0.050 

-0.62206 

-0.60598 

10 6 

0.100 

-1.03724 

-1.02373 

10 6 

0.150 

-1.31348 

-1.30265 

10® 

0.175 

-1.40878 

-1.39931 

10 6 

0.200 

-1.47881 

-1.47070 
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Fig. 5. Combining gain versus p max with L = 5000 for best-case Jfe. 
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a substantial SNR gain in this best-case scenario, particularly for the correlated noises algorithm at larger 
values of p mox - We also note that the value of p ma x cannot exceed 0.2 in order to preserve the positive 
definiteness of R B . 

Figures 6 through 8 repeat Figs. 3 through 5, respectively, using the worst-case scenario R B matrix in- 
stead, while keeping the other parameters unchanged. These examples show that this worst-case scenario 
can result in an SNR loss, and the loss nominally increases with increasing p max . Figure 6 again shows 
that convergence of the combining gain for the correlated noises algorithm to within 0.1 dB of its limiting 
value occurs at about L = 1000 samples. The small sample performance of the correlated algorithm 
again lags that of the uncorrelated algorithm. Similarly to Fig. 4, Fig. 7 shows that the performance 
superiority of the correlated noises algorithm increases with increasing (3. However, as /3 approaches one, 
the SNR loss of the correlated noises algorithm turns around and starts to decrease with increasing p max . 
At large values of (3 close to one, the correlated noises algorithm can in fact have an SNR gain even 
when the uncorrelated noises algorithm still has an SNR loss. Finally, Fig. 8 shows that the SNR loss 
of the correlated noises algorithm in this unfavorable situation deteriorates much slower than that of the 
uncorrelated noises algorithm as p ma x increases. 



NUMBER OF SAMPLES L 

Fig. 6. Combining gain versus L with p = 0.7 for worst-case Rg- 


VI. Conclusion 

The correlated noises signal-combining scheme developed in this article has been shown to be asymp- 
totically optimal as the number of residual carrier-signal samples used in the estimates of the optimum 
weight coefficients increases. The numerical examples considered here show that convergence of the com- 
bining algorithm’s SNR performance to the optimum achievable SNR performance level to within 0.1 dB 
occurs at about L = 1000 signal samples. Hence, real-time operation is possible, although the inversion 
of a complex-valued A-dimensional matrix is required at every update of the weight coefficient estimates. 
The numerical examples considered here consistently demonstrate a significant performance superiority 
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FRACTION OF TOTAL POWER 


Fig. 7. Combining gain versus P with L = 5000 for worst-case Bb- 
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MAXIMUM CORRELATION COEFFICIENT 
Fig. 8. Combining gain versus p max with L = 5000 for worst-case Rg. 


47 





of the correlated noises combining scheme over the uncorrelated noises signal-combining algorithm. It 
appears to be the combining system of choice in the presence of significant correlations between the noise 
processes in different array feed received signals. A degree of caution should be exercised in extrapolating 
the expected amount of SNR gain or loss from the numerical examples considered in Section V. The best- 
case and worst-case performances in these examples should not be viewed as being typical. Moreover, 
these two cases also should not be viewed as being best- and worst-case performances in general. 
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Appendix A 

Combiner Output Signa! Variance 


An explicit expression for Var [s c (i B )] is derived in this appendix. For convenience in the following 
analysis, let us use W M I as the shortened notation for W m l Ua), X_ml f° r 2Lml(}a), and R A ML for 
Ba,ml(^a)- Using Eqs. (26) and (22), we can write 


Var[s c (i B )] =Var [w T ML VeP a ^ s ' 


=e [(wl L r) (k t ml v)] - | W%V 


=Z f E [W ML W T ML 


V - 1 wj>v | 2 


ssV^E 


(L-K- l) 2 --i * *-1 

cos 2 # -—aml2Lml*mlBa,ml 


v-\wlv\ 2 


(A-l) 


Since R a ,ml independent of X_ ML , we can write 
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• - -i - - 1 ~ -i 

Ba,ml2Lml2LmlBa,ml 


' =E [E [r a i 


ml^ml£-mlBa,ml I Ba,ml 




-E R Ay ML E I B A ,ML Ba,ML 


— E Sa.mlE [A ML A ML j R Ay ML 


(A-2) 


Moreover, since X ml Ua) is a complex Gaussian random vector with mean X_ and covariance matrix 
(1 /L)R a , we have 


E [x ML xi L ]=^R A + XX^ 


(A-3) 


Therefore, it follows from Eqs. (A-l), (A-2), and (A-3) that 


Var[s c (i B )] = 


{L-K- 1 )' 
L 2 77 2 cos 2 <5 


U f E 


Ra.ml (jBa+XX^R 


-i 

■A, ML 


V-\W%V\ 2 


(L-K-iy 

r] 2 cos 2 6 


U f E 


(jBa+ 2UL^ (LR AML y^V- | WlV | 2 (A-4) 
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Let Tr(A) denote the trace of a square matrix A and let I_k denote the K x K identity matrix. Recall 
that LR a ml is a K x K CW(R A , L — 1) distributed random matrix. It has been shown in [8] that if A 
is a K x K CW(R a , L — 1) random matrix and C_ is any constant K x K matrix, then 


E [- 1 (L — K)(L - K -2) 


Ea'CRA 1 + 


(L - K)(L -K- 1 )(L - K - 2) 


Tr^ 1 )^ 1 


(A-5) 


for L > K. So, using Eq. (A-5) in Eq. (A-4) gets 

V ” MW1 -I L-W-K-] " *> - *1 + *?**&) 

+ Tr (±l K + XX' Ea 1 ) Sa } Y. - | W£v | 2 

Ui [i -) +*&**&£ 

+ (L-K)iAl~-L^A ,Tr (r& + “’«') ~ ' I s (a-6) 


Since A = V. cos 6, and R s = r]R A , it follows from Eq. (29) that 


^ 2 - 6 (y}R A l X2C}R A l v) = I V'R-^V I 2 = I WZ,V I 2 = 7 2 MAX (A-7) 


Since R. A l is Hermitian, Tr( AA t R 1 1 ) = XJ R A X_- So we have, by using Eq. (A-7), 




?7 2 cos 2 (5 


-r+X^RT/X 


V'Ea'V 


K 

Lrf 2 cos 2 8 


(V-'Ea'y) +v ] Eb 1 Yv}Eb 1 v 


X 2 

f7 max + 7max 


Lrjcos 2 8 

Finally, using Eqs. (A-7) and (A-8) in Eq. (A-6) yields the following expression for Var[s c (i B )]: 


(A-8) 


Var [s c (i B )] 


(L - K - 1 )(L - 1) / \ 


+ 


(L - A)(L - K - 2) \Lr]cos 2 8 ) L — K — 2 


„,2 

A MAX 


(A-9) 


Note that Var[s c (i B )] approaches zero in the limit as L tends to infinity. This is because the estimate 
W_ ml {ia) converges with probability one to the optimum weight vector W Q given by Eq. (23) as L tends 
to infinity. 
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Appendix B 

Combiner Output Noise Variance 


Consider next an explicit expression for Var[n c (is)]. We shall also employ here the shortened notations 
W ML' X Ml • ns in Appendix A. Moreover, we shall use the shortened notation n for n(*s). Since 

W ML is independent of n, a derivation similar to that establishing Eq. (A-2) can be used along with 
Eq. (22) to get 


Var[n c (* B )] =E [wl L nrSW ML 


=E [W T ML R B W* ML ] 


-J2 ~ TJCOS 2 6^ aRa,MlX-ML (B-l) 


Since R a ,ml ls independent of X.ml> using the same approach on the expected value in Eq. (B-l) gets 




K' ml E 


{^B^aml) —a {^Ra,ml) 


%-ML 


(B-2) 


Next, using the property of Eq. (A-5) for the KxK CW(R A , L — 1) distributed random matrix LR aml 
in Eq. (B-2) results in the following expression: 


Var [n c (i B )\ 


(L-K- l ) ra r*t [(L-1)-K\E?+KEa* 

~ 77cos 2 <5 b [~ ML ( L - K)(L -K-2) ~ ML , 


(L-K-1){L-1) 
>cos 2 <5(L - K){L - K - 2) 


E 


[^mlBTaILml 


It follows from Eq. (A-3) that 

E |^A ML; R a 1 A ml | =E jTr [k a X.ml^Lml}^ 
=Tr (r^E [KmlKml]) 

=j+x 1 Ra 1 2L 


(B-3) 


(B-4) 
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Finally, using Eq. (B-4) in Eqs. (B-3) and (13) yields 


Var [n c (i B )\ 


1 + 


LK - K 2 - K + 1 


(L - K)(L - K - 2) J 


7 MAX + 


(L-1)(L-K-1) 


K 


(L - K)(L - K - 2)t ? cos 2 «5 L 


(B-5) 


Appendix C 

Mean-Square Convergence of Estimated Weights 


We employ the shortened notations W_ ML , ILmli an< ^ Qa,ml> similar to the usage in the previous 
appendices. Since 


W 


ML 


W r 


= Tr 


(Kml 


Ko) (] W ML 


w, 


»)1 


and 


E 




= E 


K ML W 


* 1 

ML 


- Wn Wl 


we need only establish that E \W ML w\ iL ] — > W 0 W p as L — > oo to prove the mean-square convergence 
of W.ml to Wo- Using Eqs. (17) and (18), we can write 


E 



(L - 1 - K) 2 
L 2 r] 2 cos 2 6 



-i . * -t 
X-mlZml 



It then follows from using Eqs. (A-2) and (A-3) in Eq. (C-l) that 
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.mlW-ml 


{L-l- K)‘ 
; 2 6 


E 


Tj COS 


(lR^^ (^R* a + X*X T ^ [lR\ m ^ 1 


(C-l) 


(C-2) 


So, by using the property of Eq. (A-5) for the CW(R A ,L - 1) distributed random matrix LR aml in 
Eq. (C-2), we can write 


E 


\w ML wl L 


L-l-K 

T) 2 COS 2 8 


(L-l-K) 


\ mr 1 + mr l 2rx T (e^)- 1 ] + [f + * T (Ea)- 1 ^] mr 


(L - K)(L - K - 2) 


Therefore, E [W ML w} ML \ — > (RY)~ l X* ^(R^Y 1 lY cos 2 8 = W 0 W} 0 as L — > oo, thereby establish- 
ing the mean-square convergence of W MI to W 0 . 
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Channel Capacity of an Array System for Gaussian 
Channels With Applications to Combining and 

Noise Cancellation 
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A closed-form expression for the capacity of an array of correlated Gaussian 
channels is derived. It is shown that when signal and noise are independent, the 
array of observables can be replaced with a single observable without diminishing the 
capacity of the array channel. Examples are provided to illustrate the dependence 
of channel capacity on noise correlation for two- and three-channel arrays. 


I. Introduction 

In this article, we formulate the framework to evaluate the channel capacity of an array system. We 
define the channel capacity of an array channel as the maximum of the mutual information between a 
singl e input source and an array of n output observables over all distributions on the input that satisfy 
certain constraints (e.g., power, bandwidth, discrete or continuous, etc.). We derive a closed-form general 
formula to the channel capacity of an array of n Gaussian channels. This formula applies to correlated and 
uncorrelated noise, as long as the Gaussian assumption holds and the second-order statistics (covariance) 
of the signal and noise sources can be characterized. This formula allows one to find the channel capacity 
of various array constellations and orientations in the presence of correlated and uncorrelated noise. Some 
of the interesting results that we observed are as follows: 


(1) When the noise sources are uncorrelated, the array channel capacity is equivalent to the 
channel capacity of a single Gaussian channel with a signal-to-noise ratio (SNR) equal 
to the sum of the SNRs of the individual channels. 

(2) The array channel capacity is lower bounded by the channel capacity of the channel with 
the highest SNR. 


II. Problem Formulation 

We consider transmitting a single complex source through n channels, as illustrated in Fig. 1. Let 
a x ,ij = EXiXf,crz,ij = EZiZ ■?, and o Y ,ij = EYjY J for 1 < i, j < n. Notice that o 2 xl , o X2 , • • • , o X n are 
the signal powers as seen by the receivers; cr| v cr| 2 , • • • , o\ n are the noise variances; and Oy x , Oy 2 , • • • , 
Oy n are the channel output variances. The array channel capacity is given by 
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CHANNEL 



C array — m&X I{Y \ , Y 2 , • • • , Y n , X ) 

p(x) 

= max{ff(Fi, r 2 , ■■■,Y n )-H(Y 1 ,Y 2 ,---,Y n \X)} 

p(x) 

= max{H(Y u Y 2 , ■ ■ ■ , Y n ) - F(Z 1; Z 2 , • • • , Z„)} (1) 

p(x) 

where we further assume a power constraint on X, and where the Xi are obtained from X by a deter- 
ministic operation on each i. The interpretation of this model is that signals in the various channels 
can have different magnitudes and phases, but that the differential delays between the waveforms are 
negligible, having been removed by delay compensation. Thus, the signal path differences between the 
various channels either are negligibly small, as in the case of array-feed or phased-array applications, or 
have been properly equalized, as in the case of antenna arrays. 


HI. Capacity of an Array of Gaussian Channels 

In the following derivation, we will use some results from Cover and Thomas [1]. Since Z T = 
(Zi, Z 2 , • • • , Z n ) is a complex Gaussian random vector, H(Z) is given by [1, Theorem 9.4T] 

H(Z) — - log 2 (27re) n |©z| bits/channel use (2) 

where ©z is the covariance matrix of Z, and |©z| is its determinant. Prom Theorem 9.6.5 of [1] and 
under the assumption that the input source is power constrained to P, the input source X that maximizes 
H(Yi,Y 2 , - • ■ ,Y n ) has a Gaussian distribution. The maximum mutual information and, therefore, the 
array channel capacity are achieved with a Gaussian source, and the output observables Yi,Y 2 , • • • ,Y n 
are correlated complex Gaussian variables. Using Theorem 9.4.1 of [1], H(Yi,Y 2 , ■ ■ ■ , Y n ) is given by 
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( 3 ) 


H(Y u Y 2 , ■ ■ ■ ,Y n ) = ilog 2 (27re) n |©y | bits/channel use 


where ©y is the covariance matrix of Yj, Y 2 , ■ ■ ■ , Y n . The array channel capacity is 


y array 


1 |0y 1 

— log 2 — — - bits/channel use 

2 |©z| 


( 4 ) 


This formula can also be found in [1, Eq. (10.98)]. Notice that this formula makes no assumption on the 
pairwise correlation between the signal Xi and the noise Zj. 

Next, we consider the problem from the viewpoint of communications theory and assume that the 
additive complex Gaussian noise Z_ is independent of the signal X. For the Gaussian channel, we let 
@x = S gt, where S T = (S u S 2 , ■ ■ • , S n ) is a deterministic vector with Y^i=i |S'i| 2 = P, and f is defined 
as the conjugate transpose of a vector (that is, = S* T ). Now the covariance matrix of the observables 
can be expressed as 


©y = Q x + © z = S + Q z 


( 5 ) 


since, in this case, ESiZj = 0 for 1 < i, j < n. With the above notation, Eq. (4) can be evaluated as 


C array — 


|5 5t + e z | 

lezl 


bits/channel use 


( 6 ) 


emphasizing the independent signal and noise components of the observable covariance matrix. Using a 
corollary to Theorem A.3.2 in [3], the ratio of determinants in Eq. (6) can be written as a quadratic form 
of the inverse noise covariance matrix and the signal vector as 


| S5t +© z | 
\*z\ 


= i +jgte£ 1 £ 


( 7 ) 


and, hence, the array capacity may be equivalently expressed as 

Carray = ~ log 2 (l + /^©y 1 /?) bits/channel use 


(8) 


While Eqs. (6) and (8) are mathematically equivalent, Eq. (8) is particularly important for the following 
reasons: First, it provides useful insights into the behavior of array capacity and, second, it suggests a 
simple receiver structure for processing the array observables. 


IV. Receiver Structure for an Array of Gaussian Channels 

Let w T = (wi,w 2 , ■ ■ ■ ,w n ) be a complex weight vector. It follows that the SNR of the scalar output 
v = w T Y (Y t = (Yi, Y 2 , • • • , Y n )) is given by 
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SNR = 


\w T S\ 2 

E(\w T Z\ 2 ) 


w T SS^w* 

wQzw* 


This expression holds for any weight vector. As shown in the Appendix, the weight vector w that 
maximizes the SNR is given by 


y^opt — a (®z 1 ^)* (10) 

where a is an arbitrary complex constant. Substituting Eq. (10) into Eq. (9) yields the optimal SNR 

SNR max = ste^S ( 11 ) 

which is seen to be the same as the quadratic form in Eq. (8). Thus, the quadratic form in Eq. (8) is 
equivalent to the maximum of the SNR obtained from an optimally weighted sum of the array observables. 
The array receiver structure implied by this observation is shown in Fig. 2. 



Fig. 2. Receiver structure derived from Eq. (8). 

It is well known that a Gaussian source achieves capacity for a Gaussian channel [1]. Since the output 
of the receiver in Fig. 2 is a weighted sum of the Xi plus Gaussian noise, it is a Gaussian random variable 
for any value of the source X : hence, the output is a Gaussian channel. Since the array capacity in Eq. (8) 
depends only on the maximum SNR of the output variable v, it follows that the capacity of the scalar 
channel of Fig. 2 equals the capacity of the array channel of Fig. 1. This is an important observation since 
it enables the conversion of an n-dimensional vector observable to a single-dimensional scalar observable 
without decreasing the capacity of the channel. 

Writing Eq. (8) in terms of Eq. (11) emphasizes the point that the channel capacity of the array 
depends only on the maximum SNR that can be achieved by a weighted sum of the array observables: 

Carroty = 2 ^°§2(1 T ^^moi) (12) 

It follows that the maximum of the mutual information between v and X can also be stated directly as 
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max{H{v) - H(v \X)} = C v = C array (13) 

p{x) 

which is simply the capacity C v of the scalar channel. 

When the components of the noise vector Z_ are uncorrelated, both the covariance matrix G z and 
its inverse Q^ 1 become diagonal matrices, with the diagonal element of the inverse matrix equal to the 
inverse of the corresponding diagonal element of the covariance matrix. With the ith diagonal element of 
Qz equal to g z % the maximized SNR becomes 


SNRma, = 5t Q^S = ^ ^ (14) 

i = 1 ^Z,i 

where the right-hand side of Eq. (14) is recognized as the sum of the individual channel SNRs. The array 
capacity follows directly as 

1 / " I 0.12 

Carray = X '°g2 1 + X/ 

1 \ i= i a z,i 

Thus, the capacity of an array of Gaussian channels with uncorrelated noise is equivalent to that of a 
single Gaussian channel, with an SNR equal to the sum of the individual channel SNRs. 

For the correlated noise case, the observation of the noise in any channel provides useful information 
about the noise in every other channel. This concept is called “noise cancellation” in the adaptive signal- 
processing literature. 


bits/channel use 


(15) 


V. Examples 

Some examples of simple array channels that allow closed-form analytic solutions and provide insights 
into the problem are considered. 

A. Two-Channel Array 

Consider a two-channel array with signal vector S_= {Si, S^) and noise covariance matrix 


©z = 


<7 


2 

Z,1 


P*GZ,1&Z,2 


PVZ,l°Z,2 

4,2 J 


(16) 


where EZiZ J = paz,i&z, 2 , E\Zi\ 2 — and F|Z 2 1 2 = <j 2 Z 2 - The determinant of the noise covariance 
matrix is |©z| = o'i,i fJ i, 2 ( 1 “ |p| 2 ); its inverse is given by 


0^ = 


\P I 2 


"z, i 
-P* 

L Gz ,lCZ ,2 


GZ,\Gz,2 

1 

4,2 J 


(17) 


and the resulting array channel capacity is 
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(18) 


C, 


array 


; l°g2 1 + 


|5 X | 2 , |S 2 | 2 2 Re{pSlS 2 ) 


i - IpI 2 


<7 


+ 


Z,1 


a 


Z, 2 


az,i<?z,2 


As the magnitude of the correlation coefficient approaches 1, the array capacity grows without bound 
except for the special case when \S\\/az,i = and the phase of p cancels the phase of 

This corresponds to the singular detection case in communications theory where the signal is recovered 
without error in the absence of noise. If the noise is uncorrelated (p = 0), the array capacity depends 
only on the sum of the channel SNRs, as noted above. If the noise is correlated but the signal is 0 in one 
of the channels (for example, |5 2 | =0), the array capacity reduces to 

+ |S 2 |=0 (19) 


This “noise-cancellation” result is independent of the phase of the correlation coefficient, implying that 
the complex noise samples in the two channels need not be phase aligned for perfect cancellation — in 
effect, the optimum weights defined in Eq. (10) ensure that the noise phases are properly aligned. The 
array capacity for the two-channel noise-cancellation case corresponding to Eq. (19) is shown in Fig. 3 
for several SNRs. 


<D 

C 

C 

-C 

.o 

C/5 


CL 

< 

o 


< 

X 

o 



Fig. 3. Array channel for the two-channel case, ISgl = 0. 


B. Three-Channel Examples 

Next, we consider some three-channel examples, but only for the case of real signals and noise (in 
other words, here we ignore phase effects). We consider the triangular and the linear constellations as 
shown in Fig. 4, each having an array of three elements. We assume that the correlation coefficient of Zi 
and Zj, which is denoted by py, is geometrically related to the distance dy of array elements i and j as 
follows: 


Pa = P ldii{ (20) 

where max{— 1, a} < p < min{l, 6}, and a and b define the range [a, 6] C [—1,1] such that the covariance 
matrix ©z is legal, that is, ©z is non-negative definite and |©z| > 0. Also, 
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Fig. 4. Array configurations of three elements: (a) triangular configuration and 
(b) linear configuration. 


EZ^Z j — pij& (^ 1 ) 

The covariance matrices Q z and 6% of the triangular and linear constellations, therefore, are given by 


and 




"z,i 

perzperzp 

POz,lO-z,'i 


perzperzp 

perzperzp 


per z per z p 
per zperzp 

4,3 . 


e 2 


Z ~ 


per z per z p 

4,2 

Ip^erzperzp per z per z pi 


u zp 

perzperzp 


P erzperzp 
perzperzp 
4,3 . 


( 22 ) 


(23) 


By substituting Q z and 0| into Eq. (8), we can evaluate the array channel capacities of the above array 
constellations as a function of p. We consider two cases: channels having different SNRs and channels 
having the same SNR. 

1. Channels Having Different SNRs. Let er 2 xi = l,cr X2 = 1, and er X3 = 9; let o zl = 
0.16, a\ 2 = 0.49, and cr| 3 = 1.44. The array channel capacities of the two constellations as a function of 
p are given in Fig. 5. We observe that the array channel capacities approach infinity when p approaches 
—0.5 and 1.0 for the triangular constellation (when |©^| = 0) and when p approaches —1.0 and 1.0 for 
the linear constellation (when |0|| = 0). The channel capacities of both constellations are the same 
(3.2 bits/channel use) at p = 0, and they are lower bounded at 2.75 bits/channel use for this example. 
Thus, the maximum degradation due to noise correlation is only 0.45 bits/channel use. 

2. Channels Having the Same SNR. Let cr x l = 1.0, er 2 X 2 = 1.0, and er X 3 = 0.25; let a z l = 
0.16, cr| 2 = 0.16, and <r| 3 = 0.04. The array channel capacities of the two constellations as a function 
of p are given in Fig. 6. Again we observe that the array channel capacities approach infinity when p 
approaches —0.5 for the triangular constellation and when p approaches —1.0 for the linear constellation. 
However, both channel capacities approach 1.43 bits/channel use for the real signal and noise case, the 
channel capacity of a single channel, when p approaches 1.0. This is apparent from the fact that the 
receiver actually sees three scaled copies of the received signal plus noise, and this is equivalent to looking 
at one channel alone. 

The observations described in the two examples conform to our intuition, and the general formula 
given in Eq. (8) predicts all these observations. 
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Fig. 5. Array capacity of Gaussian channels: 
different SNRs. 



Fig. 6. Array capacity of Gaussian channels: 
same SNR. 


VI. Summary 

The capacity of an array of n Gaussian channels has been derived. The array channel was modeled 
as n observations of a single source, in the presence of additive Gaussian noise. It was shown that an 
optimally weighted sum of the array outputs achieves the same channel capacity as the array channel. 
Several examples of two- and three-channel arrays were discussed, and graphs of channel capacities were 
provided to illustrate array capacity as a function of noise correlation as well as display examples of 
singular behavior. 
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Appendix 

Derivation of Optimum Combining Weights 


The following is a simple derivation of the optimum combining weights that maximize the SNR. Other 
derivations can be found in the adaptive signal processing literature [4,5] and in this issue [6], 

We start with a derivation for the uncorrelated noise case (diagonal covariance matrix). It is shown in 
[2], using the Cauchy-Schwarz inequality, that for the uncorrelated case the optimum combining weight 
vector JJ_ is proportional to 


u = (eD 1 vy 


where V_ is a signal vector, 0£» is a diagonal matrix with components a\ i , and O^, 1 is its inverse. When 
the weight vector U_ is applied, the combined SNR is maximized, achieving its upper bound, 


SNR max = v^e- l v 


Next, consider a correlated Gaussian noise vector with covariance matrix © z , and let D be a unitary 
matrix that diagonalizes © z . With the conjugate transpose of D, D~ 1 its inverse, and £>t = D~ l 
(unitary), we can write the diagonal covariance in terms of © z and D as 


©d = dq z d t 


Thus, D rotates vectors from the uncorrelated into the correlated reference frame, without changing their 
length, while its inverse rotates them in the opposite direction. Let S_ = D~ l V_ be a signal vector, and let 
W = D~ l U_ be a weight vector in the correlated frame. The optimum weights can be written in terms of 
© z and S as 
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u = (e-^vy = (D*e* z D^)~ 1 v* 


= D^sy 


This is the optimal weight vector in the uncorrelated frame, in terms of © z x and S_. Applying the inverse 
rotation operator £> _1 , we obtain the optimum weight vector in the correlated frame as 

w = d-'u = zr ^(e^s)* = (©^s)* 
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In this article, we present two versions of a simplified maximum a posteriori 
decoding algorithm. The algorithms work in a sliding window form, like the Viterbi 
algorithm, and can thus be used to decode continuously transmitted sequences 
obtained by parallel concatenated codes, without requiring code trellis termination. 
A heuristic explanation is also given of how to embed the maximum a posteriori 
algorithms into the iterative decoding of parallel concatenated codes (turbo codes). 
The performances of the two algorithms are compared on the basis of a powerful 
rate 1/3 parallel concatenated code. Basic circuits to implement the simplified a 
posteriori decoding algorithm using lookup tables, and two fmther approximations 
(linear and threshold), with a very small penalty, to eliminate the need for lookup 
tables are proposed. 


I. Introduction and Motivations 

The broad framework of this analysis encompasses digital transmission systems where the received 
signal is a sequence of wave forms whose correlation extends well beyond T, the signaling period. There 
can be many reasons for this correlation, such as coding, intersymbol interference, or correlated fading. It 
is well known [1] that the optimum receiver in such situations cannot perform its decisions on a symbol- 
by-symbol basis, so that deciding on a particular information symbol Uk involves processing a portion of 
the received signal Td seconds long, with Td > T. The decision rule can be either optimum with respect 
to a sequence of symbols, u% = (n*,, Uk+i, ■■■ , Uk+n-i), or with respect to the individual symbol, Uk- 

The most widely applied algorithm for the first kind of decision rule is the Viterbi algorithm. In its 
optimum formulation, it would require waiting for decisions until the whole sequence has been received. 
In practical implementations, this drawback is overcome by anticipating decisions (single or in batches) 
on a regular basis with a fixed delay, D. A choice of D five to six times the memory of the received data 
is widely recognized as a good compromise between performance, complexity, and decision delay. 

Optimum symbol decision algorithms must base their decisions on the maximum a posteriori (MAP) 
probability. They have been known since the early seventies [2,3], although much less popular than the 
Viterbi algorithm and almost never applied in practical systems. There is a very good reason for this 
neglect in that they yield performance in terms of symbol error probability only slightly superior to 
the Viterbi algorithm, yet they present a much higher complexity. Only recently, the interest in these 
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algorithms has seen a revival in connection with the problem of decoding concatenated coding schemes. 
Concatenated coding schemes (a class in which we include product codes, multilevel codes, generalized 
concatenated codes, and serial and parallel concatenated codes) were first proposed by Forney [4] as a 
means of achieving large coding gains by combining two or more relatively simple “constituent” codes. 
The resulting concatenated coding scheme is a powerful code endowed with a structure that permits an 
easy decoding, like “stage decoding” [5] or “iterated stage decoding” [6]. 

To work properly, all these decoding algorithms cannot limit themselves to passing the symbols decoded 
by the inner decoder to the outer decoder. They need to exchange some kind of soft information. Actually, 
as proved by Forney [4], the optimum output of the inner decoder should be in the form of the sequence 
of the probability distributions over the inner code alphabet conditioned on the received signal, the a 
posteriori probability (APP) distribution. There have been several attempts to achieve, or at least to 
approach, this goal. Some of them are based on modifications of the Viterbi algorithm so as to obtain, at 
the decoder output, in addition to the “hard” -decoded symbols, some reliability information. This has led 
to the concept of “augmented-output,” or the list-decoding Viterbi algorithm [7], and to the soft-output 
Viterbi algorithm (SOVA) [8]. These solutions are clearly suboptimal, as they are unable to supply the 
required APP. A different approach consisted in revisiting the original symbol MAP decoding algorithms 
[2,3] with the aim of simplifying them to a form suitable for implementation [9-12]. 

In this article, we are interested in soft-decoding algorithms as the main building block of iterative stage 
decoding of parallel concatenated codes. This has become a “hot” topic for research after the successful 
proposal of the so-called turbo codes [6]. They are (see Fig. 1) parallel concatenated convolutional codes 
(PCCC) whose encoder is formed by two (or more) constituent systematic encoders joined through an 
interleaver. The input information bits feed the first encoder and, after having been interleaved by the 
interleaver, enter the second encoder. The codeword of the parallel concatenated code consists of the 
input bits to the first encoder followed by the parity check bits of both encoders. Generalizations to more 
than one interleaver are possible and fruitful [13]. 


RATE 1/3 PCCC ^ 



■J 


Fig. 1. Parallel concatenated convolutional code. 

The suboptimal iterative decoder is modular and consists of a number of equal component blocks 
formed by concatenating soft decoders of the constituent codes (CC) separated by the interleavers used 
at the encoder side. By increasing the number of decoding modules and, thus, the number of decoding 
iterations, bit-error probabilities as low as 10 -5 at E b /No = 0.0 dB for rate 1/4 PCCC have been shown 
by simulation [13]. A version of turbo codes employing two eight-state convolutional codes as constituent 
codes, an interleaver of 32 x 32 bits, and an iterative decoder performing two and one-half iterations 
with a complexity of the order of five times the maximum-likelihood (ML) Viterbi decoding of each 
constituent code is presently available on a chip yielding a measured bit-error probability of 0.9 x 10 -6 
at E b /N 0 = 3 dB [14], 
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In recent articles [15,17], upper bounds to the' ML bit-error probability of PCCCs have been proposed. 
As a by-product, it has been shown by simulation that iterative decoding can approach quite closely the 
ML performance. The iterative decoding algorithm was a simplification of the algorithm proposed in [3], 
whose regular steps and limited complexity seem quite suitable to very large-scale integration (VLSI) 
implementation. Simplified versions of the algorithm [3] have been proposed and analyzed in [12] in the 
context of a block decoding strategy that requires trellis termination after each block of bits. Similar 
simplification also was used in [16] for hardware implememtation of the MAP algorithm. 

In this article, we will describe two versions of a simplified MAP decoding algorithm that can be used 
as building blocks of the iterative decoder to decode PCCCs. A distinctive feature of the algorithms is 
that they work in a “sliding window” form, like the Viterbi algorithm, and thus can be used to decode 
“continuously transmitted” PCCCs, without requiring trellis termination and a block-equivalent structure 
of the code. The simplest among the two algorithms will be compared with the optimum block-decoding 
algorithm proposed in [3]. The comparison will be given in terms of bit-error probability when the 
algorithms are embedded into iterative decoding schemes for PCCCs. We will choose, for comparison, 
a very powerful PCCC scheme suitable for deep-space applications [18-20] and, thus, working at a very 
low signal-to-noise ratio. 


II. System Context and Notations 

As previously outlined, our final aim is to find suitable soft-output decoding algorithms for iterated 
staged decoding of parallel concatenated codes employed in a continuous transmission. The core of such 
algorithms is a procedure to derive the sequence of probability distributions over the information symbols’ 
alphabet based on the received signal and constrained on the code structure. Thus, we will start by this 
procedure and only later will we extend the description to the more general setting. 

Readers acquainted with the literature on soft-output decoding algorithms know that one burden in 
understanding and comparing the different algorithms is the spread and, sometimes, mess of notations 
involved. For this reason, we will carefully define the system and notations and then stick consistently to 
them for the description of all algorithms. 

For the first part of the article, we will refer to the system of Fig. 2. The information sequence u, 
composed of symbols drawn from an alphabet U — {u\, • • • , u/} and emitted by the source, enter an 
encoder that generates code sequences c. Both source and code sequences are defined over a time index 
set K (a finite or infinite set of integers). Denoting the code alphabet C = {ci, • • • , cm}, the code C can 
be written as a subset of the Cartesian product of C by itself K times, i.e., 

CCC K 


The code symbols Ck (the index k will always refer to time throughout the article) enter the modulator, 
which performs a one-to-one mapping of them with its signals, or channel input symbols Xk, belonging 
to the set X = {xi, • • • , xm} 1 

The channel symbols Xk are transmitted over a stationary memoryless channel with output symbols yk- 
The channel is characterized by the transitions probability distribution (discrete or continuous, according 
to the channel model) P(y\x). The channel output sequence is fed to the symbol-by-symbol soft-output 
demodulator, which produces a sequence of probability distributions 7 fc(c) over C conditioned on the 
received signal, according to the memoryless transformation 


1 For simplicity of notation, we have assumed that the cardinality of the modulator equals that of the code alphabet. In 
general, each coded symbol can be mapped in more than one channel symbol, as in the case of multilevel codes or trellis 
codes with parallel transitions. The extension is straightforward. 
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Fig. 2. The transmission system. 


7 k(c) = P(x k = x(c),y k ) = P(y k \x k = x{c))P k (c) = 7 fe ( x ) (1) 

where we have assumed to know the sequence of the a priori probability distributions of the channel input 
symbols ( P k (x ) : k 6 K ) and made use of the one-to-one mapping C —> X. 

The sequence of probability distributions 7 k (c) obtained by the modulator on a symbol- by-symbol 
basis is then supplied to the soft-output symbol decoder, which processes the distributions in order to 
obtain the probability distributions Pjt(u|y). They are defined as 


Pk(u\y) = P(u k = u|y) (2) 

The probability distributions P k (u\y) are referred to in the literature as symbol-by-symbol a posteriori 
probabilities (APP) and represent the optimum symbol-by-symbol soft output. 

Prom here on, we will limit ourselves to the case of time-invariant convolutional codes with N states, 
use the following notations with reference to Fig. 3, and assume that the (integer) time instant we are 
interested in is the fcth: 


(1) Si is the generic state at time k, belonging to the set S = {S'i, ■ • • , S jv} 

(2) S~(u') is one of the precursors of Si, and precisely the one defined by the information 
symbol u' emitted during the transition S~(u') — ► Si? 

(3) S^ (u) is one of the successors of Si, and precisely the one defined by the information 
symbol u emitted during the transition Si — > S? (u). 

(4) To each transition in the trellis, a signal x is associated, which depends on the state 
from which the transition originates and on the information symbol u determining that 
transition. When necessary, we will make this dependence explicit by writing x(u', Si) 
when the transition ends in Si and x(Si,u) when the transition originates from S{. 


III. The BCJR Algorithm 

In this section, we will restate in our new notations, without derivation, the algorithm described 
in [3], which is the optimum algorithm to produce the sequence of APP. We will call this algorithm the 


2 The state and the symbol u' uniquely specify the precursor S~ (u') in the case of the class of recursive convolutional 
encoders, like the ones we are interested in (when the largest degree of feedback polynomial represents the memory 
of a convolutional encoder). The extension to the case of feed-forward encoders and other nonconventional recursive 
convolutional encoders is straightforward. 
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Fig. 3. The meaning of notations. 


BCJR algorithm from the authors’ initials. 3 We consider first the original version of the algorithm, which 
applies to the case of a finite index set K = {1, ■ • • ,n} and requires the knowledge of the whole received 
sequence y = (yi, • • • , y n ) to work. In the following, the notations u, c, x, and y will refer to sequences 
n-symbols long, and the integer time variable k will assume the values 1, •••,«. As for the previous 
assumption, the encoder admits a trellis representation with N states, so that the code sequences c (and 
the corresponding transmitted signal sequences x) can be represented as paths in the trellis and uniquely 
associated with a state sequence s = (sq, • ■ • , s n ) whose first and last states, ,s 0 and s n , are assumed to 
be known by the decoder. 4 


Defining the a posteriori transition probabilities from state Si at time k as 


a k (Si,u) = P{u k = u,s fe — i = 5j|y) 


( 3 ) 


the APP P(u|y) we want to compute can be obtained as 


PfcMy) = ^ ^a k {Si,u ) 

Si 


( 4 ) 


Thus, the problem of evaluating the APP is equivalent to that of obtaining the a posteriori transition 
probabilities defined in Eq. (3). In [3], it was proven that the APP can be computed as 


<?k{Si,u) = h a a k -i(Si)ik(.x(Si,u))Pk{Si~ (u)) 


( 5 ) 


where 


3 The algorithm is usually referred to in the recent literature as the “Bahl algorithm” ; we prefer to credit all the authors: 
L. R. Bahl, J. Cocke, F. Jelinek, and J. Raviv. 

4 Lower-case s k denotes the states of a sequence at time k, whereas upper-case Si represents one particular state belonging 
to the set S. 
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« h a is such that 


]> ~2(?k(Si,u ) = 1 

Si 


• 7 k(x(Si,u)) are the joint probabilities already defined in Eq. (1), i.e., 


7 k(x) = P(yk, x k = x) = P{y k \x k = x) ■ P{x k = x) (6) 

The 7’s can be calculated from the knowledge of the a priori probabilities of the channel 
input symbols x and of the transition probabilities of the channel P(y k \x k = x). For 
each time k, there are M different values of 7 to be computed, which are then associated 
to the trellis transitions to form a sort of branch metrics. This information is supplied 
by the symbol-by-symbol soft-output demodulator. 

• a k {Si) are the probabilities of the states of the trellis at time k conditioned on the past 
received signals, namely, 


a k (Si) = P(s k = Si\ V i) (7) 

where y f denotes the sequence • • • ,y k . They can be obtained by the forward 

recursion 5 


a k (Si) = h a J2 a k-i(Si (u))7fe(x(u,5i)) (8) 

U 


with h a a constant determined through the constraint 

£<**($) = ! 

Si 


and where the recursion is initialized as 


ao {Si) = { q 


if Si — s 0 
otherwise 


(9) 


• f3 k (Si) are the probabilities of the trellis states at time k conditioned on the future 
received signals P(s k = 5j|y£ +1 ). They can be obtained by the backward recursion 


Pk(Si) = h /3 £/?* :+ i(5+(u))7 fc+ i(x(S' i ,u)) (10) 

U 


5 For feed-forward encoders and nonconventional recursive convolutional encoders like G(D) = [1, (1 + D + D 2 )/( 1 4- D)} 
in Eq. (8), the summation should be over all possible precursors S'r (u) that lead to the state Si, and x(u, Si) should be 
replaced by x(S~ (u),u). Then such modifications are also required for Eqs. (18) and (26). In Eqs. (22), (29), and (32), 
the maximum should be over all S~( u) that lead to Si. The c{u,Si ) should be replaced by c(Sr (u), u). 
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with hp a constant obtainable through the constraint 


E&( 5 «) = 1 

Si 


and where the recursion is initialized as 


/?«($) = {J 


if Si = s n 
otherwise 


( 11 ) 


We can now formulate the BCJR algorithm by the following steps: 


(1) Initialize ao and f3 n according to Eqs. (9) and (11). 

(2) As soon as each term yk of the sequence y is received, the demodulator supplies to the 
decoder the “branch metrics” 7 of Eq. (6), and the decoder computes the probabilities 

according to Eq. (8). The obtained values of ak(Si) as well as the 7 k are stored for 
all k, Si, and x. 

(3) When the entire sequence y has been received, the decoder recursively computes the 
probabilities j3k according to the recursion of Eq. (10) and uses them together with the 
stored a’s and 7’s to compute the a posteriori transition probabilities <Jk(Si, u) according 
to Eq. (5) and, finally, the APP Pk{u |y) from Eq. (4). 


A few comments on the computational complexity of the finite-sequence BCJR algorithm can be found 
in [3]. 


IV. The Sliding Window BCJR (SW-BCJR) 

As previous description made clear, the BCJR algorithm requires that the whole sequence has been 
received before starting the decoding process. In this aspect, it is similar to the Viterbi algorithm in its 
optimum version. To apply it in a PCCC, we need to subdivide the information sequence into blocks, 6 
decode them by terminating the trellises of both CCs, 7 and then decode the received sequence block by 
block. Beyond the rigidity, this solution also reduces the overall code rate. 


A more flexible decoding strategy is offered by a modification of the BCJR algorithm in which the 
decoder operates on a fixed memory span, and decisions are forced with a given delay D. We call this 
new, and suboptimal, algorithm the sliding window BCJR (SW-BCJR) algorithm. We will describe two 
versions of the sliding window BCJR algorithm that differ in the way they overcome the problem of 
initializing the backward recursion without having to wait for the entire sequence. We will describe the 
two algorithms using the previous step description suitably modified. Of the previous assumptions, we 
retain only that of the knowledge of the initial state s 0 , and thus assume the transmission of semi-infinite 
code sequences, where the time span K ranges from 1 to 00. 


6 The presence of the interleaver naturally points toward a block length equal to the interleaver length. 

7 The termination of trellises in a PCCC has been considered a hard problem by several authors. As shown in [13], it is, 
indeed, quite an easy task. 
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A. The First Version of the Sliding Window BCJR Algorithm (SW1-BCJR) 

Here are the steps: 


(1) Initialize ao according to Eq. (9). 

(2) Forward recursion at time k: Upon receiving yk, the demodulator supplies to the de- 
coder the M distinct branch metrics, and the decoder computes the probabilities otk(Si) 
according to Eqs. (6) and (8). The obtained values of ctk(Si) are stored for all Si, as well 
as the 7fc(a:). 

(3) Initialization of the backward recursion (k > D): 

Pk(Sj) = ak(Sj), VSj (12) 


(4) Backward recursion: It is performed according to Eq. (10) from time k — 1 back to time 
k-D. 

(5) The a posteriori transition probabilities at time k — D are computed according to 

Vk-D(Si,u) = h a ■ ak-D-i(Si)lk-D(x(Si,u))(3k-D(Sf(u)) (13) 


(6) The APP at time k — D is computed as 


Pfe-cHy) = y ^o’k-pjSiju) 
Si 


(14) 


For a convolutional code with parameters (k 0 ,no), number of states N, and cardinality of the code 
alphabet M = 2 n ° , the SW1-BCJR algorithm requires storage of IV x D values of a’s and M x D values 
of the probabilities 7 k{x) generated by the soft demodulator. Moreover, to update the a’s and f3’s for each 
time instant, the algorithm needs to perform M x 2 k ° multiplications and N additions of 2 fe ° numbers. 
To output the set of APP at each time instant, we need a H-times long backward recursion. Thus, the 
computational complexity requires overall 


• (D + 1 )M x 2 k ° multiplications 

• (D + 1)M additions of 2 k ° numbers each 

As a comparison, 8 the Viterbi algorithm would require, in the same situation, M x 2 k ° additions and 
M x 2 fe °-way comparisons, plus the trace-back operations, to get the decoded bits. 

B. The Second, Simplified Version of the Sliding Window BCJR Algorithm (SW2-BCJR) 

A simplification of the sliding window BCJR that significantly reduces the memory requirements 
consists of the following steps: 


8 Though, indeed, not fair, as the Viterbi algorithm does not provide the information we need. 
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(1) Initialize c*o according to Eq. (9). 

(2) Forward recursion ( k > D): If k > D, the probabilities a k -D-i(Si) are computed 
according to Eq. (8). 

(3) Initialization of the backward recursion (k > D): 

Pk(Sj) = 1 WSj (15) 


(4) Backward recursion (k > D): It is performed according to Eq. (10) from time k — 1 back 
to time k — D. 

(5) The a posteriori transition probabilities at time k — D are computed according to 

<Tk-D(Si, u) — h a ■ ak-D-l(Si)lk-D(x(Sii u))p k - D (S ? (it)) (16) 


(6) The APP at time k — D is computed as 

iVr>H y) = a k-o(Si , it) (17) 

St 

This version of the sliding window BCJR algorithm does not require storage of the N x D values of 
a’s as they are updated with a delay of D steps. As a consequence, only N values of a’s and M x D 
values of the probabilities 7 k (x) generated by the soft demodulator must be stored. The computational 
complexity is the same as the previous version of the algorithm. However, since the initialization of the 
P recursion is less accurate, a larger value of D should be set in order to obtain the same accuracy on 
the output values Pfc_/j(u|y). This observation will receive quantitative evidence in the section devoted 
to simulation results. 


V, Additive Algorithms 

A. The Log-BCJR 

The BCJR algorithm and its sliding window versions have been stated in multiplicative form. Owing 
to the monotonicity of the logarithm function, they can be converted into an additive form passing to 
the logarithms. Let us define the following logarithmic quantities: 


r fc (x) = log [7(2)] 


A k {Si) =log[a fe (Si)] 
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Bk(Si) =log[&($)] 


E k (Si,u) = log[crfe(iSj, u)] 


These definitions lead to the following A and B recursions, derived from Eqs. (8), (10), and (5): 


Ak(Si) = log 


^expfAfc-i^ (u)) + T k (x(u,Si))} 


+ H a 


Bk(Si) = log 


exp {r fc+ i ( x(Si , u )) 4- B k + 1 (Sf (u)) } 


+ H b 


s k(Si,u) =A k _ 1 (S i ) +r k (x(Si,u)) + B k (S+(u)) + 1T S 


with the following initializations: 


( 18 ) 


(19) 

( 20 ) 


A,(Si)={° 

l — oo 

Bi(5 i )={° 

1 -oo 


if Si = s 0 
otherwise 

if Si = s n 
otherwise 


B. Simplified Versions of the Log-BCJR 

The problem in the recursions defined for the log-BCJR consists of the evaluation of the logarithm of 
a sum of exponentials: 


log 


^exp{Aj} 


An accurate estimate of this expression can be obtained by extracting the term with the highest expo- 
nential, 


Am = max A* 
i 


so that 


log 


^%xp{Ai} 


= A M + log 


1 + ex P(^i - a m } 


( 21 ) 


and by computing the second term of the right-hand side (RHS) of Eq. (21) using lookup tables. Further 
simplifications and the required circuits for implementation are discussed in the Appendix. 
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However, when Am A*, the second term can be neglected. This approximation leads to the additive 
logarithmic-BCJR (AL-BCJR) algorithm: 

A k (Si) = max [A k -i(S~ ( u )) + T k (x(u, 5*))] + H A (22) 

B k (Si) =max[B fc+ i(5t( u )) +r fe+ i(a;(5 , i ,w))] + H B (23) 

Z k (Si,u) =A k .. 1 (S i ) + r k (x(Si,u)) + B k (S+(u)) + Hz (24) 


with the same initialization of the log-BCJR. 

Both versions of the SW-BCJR algorithm described in the previous section can be used, with obvious 
modifications, to transform the block log-BCJR and the AL-BCJR into their sliding window versions, 
leading to the SW-log-BCJR and the SWAL1-BCJR and SWAL2-BCJR algorithms. 


VI. Explicit Algorithms for Some Particular Cases 

In this section, we will make explicit the quantities considered in the previous algorithms’ descriptions 
by making assumptions on the code type, modulation format, and channel. 

A. Rate 1/n Binary Systematic Convolutional Encoder 

In this section, we particularize the previous equations in the case of a rate 1/n binary systematic 
encoder associated to n binary-pulse amplitude modulation (PAM) signals or binary phase shift keying 
(PSK) signals. 

The channel symbols x and the output symbols from the encoder can be represented as vectors of n 
binary components: 


C — [cj, * * * , Cnjcj £ {0, 1} 

x = [xi, • ■ • , x n )xi 6 {A, -A} 

%k ~ [a*l, ’ * * j^fen] 

Vk = [l/AUi ' ' ' j Ukn] 

where the notations have been modified to show the vector nature of the symbols. The joint probabilities 
7 fc(x), over a memoryless channel, can be split as 


n 

7fe(^) = P{ykm\^km = •^rn)B{p^km = 3-m) (25) 

m= 1 


Since in this case the encoded symbols are n-tuple of binary symbols, it is useful to redefine the input 
probabilities, 7 , in terms of the likelihood ratios: 
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A 


^ kn 


P(jjkm\%km — 

P {y km\^ km = Al) 


\ A Pip' km — A) 

km ~ P(x km = -A) 


so that, from Eq. (25), 


7 k(x) 


TT (A krn) Cm 

m-1 ^ ^ km 

771 = 1 


(A£J C ™ 

1 + A 


= ^II [ X *m-\£m] C ' n 

m= 1 


where h 7 takes into account all terms independent of x. 
The BCJR can be restated as follows: 


a k{Si) — h-yh a ctfc— i (iSj ( u )) [Afcm ■ Ajj m ] 


m=l 


&($) - M/sE/WOW) n [W)m • A 

u m=l 

n 

cr k (Si,u) = h^h a ak-i(Si) [] [A fem ■ A^ m ] Cm( “’ Si) 0 k (S+(u)) 


\A l c m(u,Si) 
* A km\ 

(26) 

" ^m{Si ,11) 


l)m ‘ A(fc + j) m 

(27) 

(3 k (S+(u)) 

(28) 


m= 1 


whereas its simplification, the AL-BCJR algorithm, becomes 


A k (Si) = max l A k _i(S i (u)) + 


E 


Cm{ u tSi) {Akm + Afc m ) 


+ #A 


(29) 


B k {Si) = max I Bfc+i^ (u)) + ^ c m (S*, u ) (A fcm + A^ m ) 

1 777—1 

n 

E k (Si,u) =A fc _i(5j) + ^ c m (S' i ,u) (A^ + A^ m ) + B k (SA( u )) (31) 

777=1 



where A stands for the logarithm of the corresponding quantity A. 

B, The Additive White Gaussian Noise Channel 

When the channel is the additive white Gaussian noise (AWGN) channel, we obtain the explicit 
expression of the log-likelihood ratios A k i as 
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A ki = log 


P{Vki\xkt = A) 

P{yki\^ki = A) 


= log 


jp - 4 > 2 ) 

VSp eXp{ -2^ (te + ^ } 


2A 


Hence, the AL-BCJR algorithm assumes the following form: 


A k (Si) — HU1X (u)) "I" ^ ' Cmiu, Si) ^ 2/fcm "1“ A fcm ^ ^ + JifyX 


nax < i 

u ^ 


(2 A 


Bk(Si) = max ^ B k+ i(Sf ( u )) + ]T c^S*, u ) ^-^-2/fem + A^ m J > + H B 


(32) 


(33) 


771=1 
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^k{^i,u) =A k —l{Si) + ^ ' Cm(Sj i ^) f "1” Afc,nJ + P k(B+(u)) 

m — 1 ' ' 


(34) 


In the examples presented in Section VIII, we will consider turbo codes with rate 1/2 component 
convolutional codes transmitted as binary PAM or binary PSK over an AWGN channel. 


VII. Iterative Decoding of Parallel Concatenated Convolutional Codes 

In this section, we will show how the MAP algorithms previously described can be embedded into 
the iterative decoding procedure of parallel concatenated codes. We will derive the iterative decoding 
algorithm through suitable approximations performed on maximum-likelihood decoding. The description 
will be based on the fairly general parallel concatenated code shown in Fig. 4, which employs three 
encoders and three interleavers (denoted by 7r in the figure). 

Let u k be the binary random variable taking values in {0, 1}, representing the sequence of information 
bits u = (ui, • • • , u n ). The optimum decision algorithm on the fcth bit u k is based on the conditional 
log-likelihood ratio L k '- 


L k = log 


P(u k = l|y) 
P(u k = 0|y) 


= log 


= log 


£u:« fc = l P (yl U ) P ( U j) , P(u k = 1) 

£u:„ fc =o P( y |u) E u p K) og P(u k = 0) 

Eu:« fc = l P (yl X ( U ))n &k P ( U j) . P(u k = 1 ) 
£u:u fc =0 P (y|x( u )) n,yfc P ( U j) ° S P ( U k = 0 ) 


where, in Eq. (35), P(uj) are the a priori probabilities. 


( 35 ) 
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Fig. 4. Parallel concatenation of three convolutional codes. 

If the rate k 0 /n 0 constituent code is not equivalent to a punctured rate 1 /n' 0 code or if turbo trellis- 
coded modulation is used, we can first use the symbol MAP algorithm as described in the previous 
sections to compute the log-likelihood ratio of a symbol u = m, 112, ■ • • , given the observation y as 


A(u) = log 


p (u|y) 

-P(°|y) 


where 0 corresponds to the all-zero symbol. Then we obtain the log-likelihood ratios of the jth bit within 
the symbol by 


Y' , e A(u) 

L ("i) = lo S ^ " ■ f; A(u) 

2-<u:Uj=0 e 

In this way, the turbo decoder operates on bits, and bit, rather than symbol, interleaving is used. 

To explain the basic decoding concept, we restrict ourselves to three codes, but extension to several 
codes is straightforward. In order to simplify the notation, consider the combination of the permuter 
(interleaver) and the constituent encoder connected to it as a block code with input u and outputs x*, 
i = 0, l,2,3(xo = u) and the corresponding received sequences as y*, i = 0,1,2, 3. The optimum bit 
decision metric on each bit is (for data with uniform a priori probabilities) 


Eu:u fc =i p (yol u ) f, (yil u ) p (y2l u ) f> (y3lu) 

Eu:« fe =o p (yol u ) p (yil u ) p (y2|u) p (y3|u) 


but, in practice, we cannot compute Eq. (36) for large n because the permutations 7r2,7r3 imply that y 2 
and y 3 are no longer simple convolutional encodings of u. Suppose that we evaluate P(yj|u), i = 0,2,3 
in Eq. (36) using Bayes’ rule and using the following approximation: 
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( 37 ) 


P(u|yi) « JJ Pi(ufc) 

k= 1 


Note that P(u|yi) is not separable in general. However, for i = 0, P(u|yo) is separable; hence, Eq. (37) 
holds with equality. So we need an algorithm that approximates a nonseparable distribution P(u|y,) = P 
with a separable distribution P( w k) = Q- The best approximation can be obtained using the 

Kullback cross-entropy minimizer, which minimizes the cross-entropy H(Q,P) — P{log(Q/P)} between 
the input P and the output Q. 

The MAP algorithm approximates a nonseparable distribution with a separable one; however it is 
not clear how good it is compared with the Kullback cross-entropy minimizer. Here we use the MAP 
algorithm for such an approximation. In the iterative decoding, as the reliability of the {uk } improves, 
intuitively one expects that the cross-entropy between the input and the output of the MAP algorithm 
will decrease, so that the approximation will improve. If such an approximation, i.e., Eq. (37), can be 
obtained, we can use it in Eq. (36) for i = 2 and i = 3 (by Bayes’ rule) to complete the algorithm. 

Define Lik by 


Pi(u k ) = 


g**k£«/fc 

1 4 . 


(38) 


where Uk € {0, 1}. To obtain {p} or, equivalently, {Lik}, we use Eqs. (37) and (38) for i = 0,2,3 (by 
Bayes’ rule) to express Eq. (36) as 

Lk — /(yi)Lo,i<2)L3> A;) + Lok + L2k + L^k (39) 


where Zofc = 2Ayofc/a 2 (for binary modulation) and 


/(yi,Lo,L 2 ,H 3 ,A:) = log 


Eu:. fc =i p (yil”) II 

Eu:„ fc =0 P (yil U ) EU e“Aio i +L 2 . J +L 3i ) 


We can use Eqs. (37) and (38) again, but this time for i = 0, 1,3, to express Eq. (36) as 


(40) 


Lk — /(y 2 ,Lo,Li,L 3 ,fe) + Lofc + p ik + L^k 


(41) 


and similarly, 


Lk — /(y3, Lo, k) + Lok + L^ + 
A solution to Eqs. (39), (41), and (42) is 

Tik =/(yi,Lo,L 2 ,L 3 ,fc) 

Lik =/(y2,i'0,Li,L 3 , k) 

Lsk =/(y3) ^Oi Ti,L 2 , k) j 


(42) 


( 43 ) 
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for k = 1, 2, • • • , n, provided that a solution to Eq. (43) does indeed exist. The final decision is then based 
on 


Lk — Lok + Lik + L/ 2 k + L 3 k (44) 

which is passed through a hard limiter with zero threshold. We attempted to solve the nonlinear equations 
in Eq. (43) for hi, L 2 , and L3 by using the iterative procedure 

^ik +1) = a^f(yi,ioA m) ,^r\k) (45) 

for k = 1, 2, • • • , n, iterating on m. Similar recursions hold for and L^\ 

We start the recursion with the initial condition = Lo- For the computation of 

/(■), we can use any MAP algorithm as described in the previous sections, with permuters (direct and 
inverse) where needed; call this the basic decoder D iy i = 1,2,3. The L^\i = 1,2,3 represent the 
extrinsic information. The signal flow graph for extrinsic information is shown in Fig. 5 [13], which is a 
fully connected graph without self-loops. Parallel, serial, or hybrid implementations can be realized based 
on the signal flow graph of Fig. 5 (in this figure yo is considered as part of yi). Based on our equations, 
each node’s output is equal to internally generated reliability L minus the sum of all inputs to that node. 
The BCJR MAP algorithm always starts and ends at the all-zero state since we always terminate the 
trellis as described in [13]. We assumed 7Ti = I identity; however, any 7Ti can be used. 



The overall decoder is composed of block decoders Di connected in parallel, as in Fig. 6 (when the 
switches are in position P), which can be implemented as a pipeline or by feedback. A serial imple- 
mentation is also shown in Fig. 6 (when the switches are in position S). Based on [13, Fig. 5], a serial 
implementation was proposed in [21]. For those applications where the systematic bits are not transmit- 
ted or for parallel concatenated trellis codes with high-level modulation, we should set Lo — 0. Even 
in the presence of systematic bits, if desired, one can set Lo = 0 and consider yo as part of yi. If the 
systematic bits are distributed among encoders, we use the same distribution for yo among the received 
observations for MAP decoders. 

At this point, further approximation for iterative decoding is possible if one term corresponding to 
a sequence u dominates other terms in the summation in the numerator and denominator of Eq. (40). 
Then the summations in Eq. (40) can be replaced by “maximum” operations with the same indices, i.e., 
replacing Y^ u -u k =i with f° r * = 0, 1. A similar approximation can be used for i 2 fc and L 3 k in 

Eq. (43). This suboptimal decoder then corresponds to an iterative decoder that uses AL-BCJR rather 
than BCJR decoders. As discussed, such approximations have been used by replacing with max in the 
log-BCJR algorithm to obtain AL-BCJR. Clearly, all versions of SW-BCJR can replace BCJR (MAP) 
decoders in Fig. 6. 

For turbo codes with only two constituent codes, Eq. (45) reduces to 
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Fig. 6. Iterative decoder structure for three parallel concatenated codes. 


X( 1 J +1) = a< m) /(yi,Lo,i4 m) ,fc) 

4r l) = 4 m) /(y2,L 0 ,L ( 1 m) ,fc) 

for k = 1, 2, • • • , n, and m — 1, 2, - • where, for each iteration, and can be optimized (simulated 
annealing) or set to 1 for simplicity. The decoding configuration for two codes is shown in Fig. 7. In this 
special case, since the paths in Fig. 7 are disjointed, the decoder structure can be reduced to a serial mode 
structure if desired. If we optimize a[ m> and a ^ 1 , our method for two codes is similar to the decoding 
method proposed in [6], which requires estimates of the variances of Lu t and L^k for each iteration in 
the presence of errors. It is interesting to note that the concept of extrinsic information introduced in 
[6] was also presented as “partial factor” in [22], However, the effectiveness of turbo codes lies in the 
use of recursive convolutional codes and random permutations. This results in time-shift-varying codes 
resembling random codes. 

In the results presented in the next section, we will use a parallel concatenated code with only two 
constituent codes. 
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DECODED BITS 


Fig. 7. Iterative decoder structure for two parallel concatenated codes. 


VIII. Simulation Results 

In this section, we will present some simulation results obtained applying the iterative decoding algo- 
rithm described in Section VII, which, in turn, uses the optimum BCJR and the suboptimal, but simpler, 
SWAL2-BCJR as embedded MAP algorithms. All simulations refer to a rate 1/3 PCCC with two equal, 
recursive convolutional constituent codes with 16 states and generator matrix 


G{D) = 


1 + D + D 3 + D 4 
1 + D 3 + D 4 


and an interleaver of length 16,384 designed according to the procedure described in [13], using an 
S-random permutation with S = 40. Each simulation run examined at least 25,000,000 bits. 

In Fig. 8, we plot the bit-error probabilities as a function of the number of iterations of the decoding 
procedure using the optimum block BCJR algorithm for various values of the signal-to-noise ratio. It can 
be seen that the decoding algorithm converges down to BER = 10“ 5 at signal-to-noise ratios of 0.2 dB 
with nine iterations. The same curves are plotted in Fig. 9 for the case of the suboptimum SWAL2-BCJR 
algorithm. In this case, 0.75 dB of signal-to-noise ratio is required for convergence to the same BER and 
with the same number of iterations. 

In Fig. 10, the bit-error probability versus the signal-to-noise ratio is plotted for a fixed number 
(5) of iterations of the decoding algorithm and for both optimum BCJR and SWAL2-BCJR MAP de- 
coding algorithms. It can be seen that the penalty incurred by the suboptimum algorithm amounts 
to about 0.5 dB. This figure is in agreement with a similar result obtained in [12], where all MAP 
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Fig. 8. Convergence of turbo coding: bit-error probability 
versus number of iterations for various E^Nq using the 
SW2-BCJR algorithm. 



Fig. 9. Convergence of turbo coding: bit-error probability 
versus number of iterations for various E^Nq using the 
SWAL2-BCJR algorithm. 


algorithms were of the block type. The penalty is completely attributable to the approximation of the 
sum of exponentials described in Section V.B. To verify this, we have used a SW2-BCJR and compared 
its results with the optimum block BCJR, obtaining the same results. 

Finally, in Figs. 11 and 12, we plot the number of iterations needed to obtain a given bit-error prob- 
ability versus the bit signal-to-noise ratio, for the two algorithms. These curves provide information on 
the delay incurred to obtain a given reliability as a function of the bit signal-to-noise ratio. 
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Fig. 10. Bit-error probability as a function of the bit signal-to- 
noise ratio using the SW2-BCJR and SWAL2-BCJR algorithms 
with five iterations. 



Et/No 


Fig. 11. Number of iterations to achieve several bit-error 
probabilities as a function of the bit signal-to-noise ratio using 
the SWAL2-BCJR algorithm. 


IX. Conclusions 

We have described two versions of a simplified maximum a posteriori decoding algorithm working in 
a sliding window form, like the Viterbi algorithm. The algorithms can be used as a building block to 
decode continuously transmitted sequences obtained by parallel concatenated codes, without requiring 
code trellis termination. A heuristic explanation of how to embed the maximum a posteriori algorithms 
into the iterative decoding of parallel concatenated codes was also presented. Finally, the performances 
of the two algorithms were compared on the basis of a powerful rate 1 /3 parallel concatenated code. 
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Fig. 12. Number ot iterations to achieve several bit-error 
probabilities as a function of the bit signal-to-noise ratio using 
the SW2-BCJR algorithm. 
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Appendix 

Circuits to Impiement the MAP Algorithm for Decoding 
Rate 1/n Component Codes of a Turbo Code 


In this appendix, we show the basic circuits required to implement a serial additive MAP algorithm 
for both block log-BCJR and SW-log-BCJR. Extension to a parallel implementation is straightforward. 
Figure A-l shows the implementation 9 of Eq. (18) for the forward recursion using a lookup table for 
evaluation of log(l + e~ x ), and subtraction of maXj{Afc(Sj)} from A k (Si) is used for normalization to 
prevent buffer overflow. 10 The circuit for maximization can be implemented simply by using a comparator 
and selector with feedback operation. Figure A-2 shows the implementation of Eq. (19) for the backward 
recursion, which is similar to Fig. A-l. A circuit for computation of log(Pfc(u|y)) from Eq. (4) using 
Eq. (20) for final computation of bit reliability is shown in Fig. A-3. In this figure, switch 1 is in position 1 
and switch 2 is open at the start of operation. The circuit accepts Efc(5j, u) for i — 1, then switch 1 moves 
to position 2 for feedback operation. The circuit performs the operations for i — 1, 2, • • • , N. When the 
circuit accepts Ejt(5j,u) for i = N, switch 1 goes to position 1 and switch 2 is closed. This operation is 
done for u = 1 and u = 0. The difference between log(P fc (l|y)) and log(P fc (0|y)) represents the reliability 
value required for turbo decoding, i.e., the value of L k in Eq. (35). 



NORMALIZED A k (Sj) 


Fig. A-1. Basic structure for forward computation in the log-BCJR MAP algorithm. 


9 For feed-forward and nonconventional recursive convolutional codes, the notations in Fig. A-l should be changed according 
to Footnotes 2 and 5. 

10 Simpler normalization can be achieved by monitoring the two most significant bits. When both of them are one, then we 
reset all the most significant bits to zero. This method increases the bit representation by an additional 2 bits. 
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NORMALIZED Bfc(S/) 


Fig. A-2. Basic structure for backward computation in the log-BCJR MAP algorithm. 


We propose two simplifications to be used for computation of log(l + e x ) without using a lookup 
table. 

Approximation 1: We used the approximation log(l + e~ x ) « —ax + b, 0 < x < b/a where b = log(2), 
and we selected a = 0.3 for the simulation. We observed about a 0.1-dB degradation compared with the 
full MAP algorithm for the code described in Section VIII. The parameter a should be optimized, and it 
may not necessarily be the same for the computation of Eq. (18), Eq. (19), and log(Pfc(u|y)) from Eq. (4) 
using Eq. (20). We call this “linear” approximation. 

Approximation 2: We take 


log(l + e x 



0 

c 


if x > T] 
if x < rj 


We selected c = log(2) and the threshold 77 = 1.0 for our simulation. We observed about a 0.2-dB 
degradation compared with the full MAP algorithm for the code described in Section VIII. This threshold 
should be optimized for a given SNR, and it may not necessarily be the same for the computation 
of Eq. (18), Eq. (19), and log(Pfc(u|y)) from Eq. (4) using Eq. (20). If we use this approximation, 
the log-BCJR algorithm can be built based on addition, comparison, and selection operations without 
requiring a lookup table, which is similar to a Viterbi algorithm implementation. We call this “threshold” 
approximation. At most, 8- to 10-bit representation suffices for all operations (see also [12] and [16]). 
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Progressive Transmission and 
Compression of Images 

A. B. Kiely 

Communications Systems and Research Section 

We describe a n image data compression strategy featuring progressive trans- 
mission. The method exploits subband coding and arithmetic coding for compres- 
sion. We analyze the Laplacian probability density, which closely approximates 
the statistics of individual subbands, to determine a strategy for ordering the com- 
pressed subband data in a way that improves rate-distortion performance. Results 
are presented for a test image. 


I. Introduction 

An image data compression system that uses progressive transmission is one that allows a user to 
reconstruct successively higher fidelity versions of an image as data are received. The goal of progressive 
transmission is thus not only efficient overall compression, but efficient compression at every step. 

If the data rate available for image transmission is unexpectedly low, or if the volume of compressed 
data exceeds expectations, the available rate will be used to its full extent, to provide nearly the highest- 
fidelity image possible given the rate constraint. Alternatively, if the available rate exceeds expectations, 
it will be possible to send higher-resolution images than originally planned. In this sense, progressive 
transmission strategies are robust with respect to the available data rate. 

In situations where the reverse channel can be used, data compression can be combined with advanced 
communications strategies to increase the volume of data returned. The use of retransmission schemes is 
one example of such a strategy [8]. Progressive transmission gives another method, because it provides 
the ability to quickly view low- or medium-resolution previews of an image, making efficient browsing 
possible. A user can decide whether to continue the transmission of the full image or proceed to the next 
image. This makes more efficient use of the channel because images of less value are not transmitted at 
full resolution. This is particularly beneficial to deep-space missions with high data volume and severe 
rate constraints. 


II. Progressive Transmission and Rate-Distortion Theory 

For a given source and distortion metric, we are interested in the trade-off between the rate (usually 
measured in average number of bits/sample) and distortion. Some of the rate-distortion functions of 
interest include 
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(1) The rate-distortion limit. This is the theoretical minimum distortion as a function of 
the average number of bits used to describe the source [2]. This describes the opti- 
mum performance obtainable, ignoring constraints such as speed and complexity, and is 
computable only if the source statistics are known. 

(2) The rate-distortion performance achievable by a particular technique. 

(3) The rate-distortion performance “progressively” achievable by a particular progressive 
transmission technique. This is the performance obtained by measuring the rate and 
distortion of the reconstructed versions of the data at each stage of the transmission. 

In theory, we can select a “target” rate-distortion limit point that represents the rate and distortion 
when all of the compressed data are transmitted. Ideally, using progressive transmission, at every point 
in the transmission we would like to be as close as possible to the rate-distortion limit. Unfortunately, it 
is not always possible to meet the target and have all of the points leading up to the target lie along the 
rate-distortion limit [4]. 

Equitz and Cover showed that it is possible to meet the rate-distortion limit progressively when 
solutions to the rate-distortion problem can be written as a Markov chain [4]. Such a source is called 
“successively refinable.” Whether or not a source is successively refinable depends not only on the source 
statistics but also on the distortion metric selected. For example, a Laplacian source is successively 
refinable with respect to mean absolute error but, in general, not with respect to other metrics [4], Even 
if a source is successively refinable, the optimal strategy for refinement may be unwieldy. Fortunately, 
even when a source is not successively refinable, the penalty for progressive transmission may be small. 

An example follows: Under mean square error (MSE) distortion, a memoryless Gaussian source with 
variance a 2 has a rate-distortion limit MSE = a 2 2~ 2R , where R is the rate in bits [2, p. 99]. In a practical 
source coding system, we might be willing to sacrifice performance in exchange for the simplicity of using, 
say, a 6-bit uniform scalar quantizer, where each quantizer interval is identified by a 6-bit codeword, and 
each codeword is compressed by arithmetic coding. For the Gaussian source, this produces a different 
rate-distortion curve parameterized by the quantizer step size. Both rate-distortion functions are shown 
in Fig. 1. 

Suppose we wish to transmit 24 independent and identically distributed Gaussian samples quantized 
using the 6-bit uniform quantizer, and we select a target rate of 4 bits/sample, which corresponds to 



0 1 2 3 4 5 


RATE, bits/sample 


Fig. 1. Some rate-distortion functions associated with a 
quantized Gaussian source. 
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MSE/a 2 « 5.56 x 10~ 3 . Figure 1 shows the progressive rate-distortion performance obtained if we 
transmit the arithmetic encoded 6-bit quantized samples sequentially. Another option is to reorder the 
codeword bits before arithmetic coding so that the 24 most-significant bits of the samples are compressed 
and transmitted first, followed by the next most significant bit from each codeword, and so on. Figure 1 
illustrates that although this bit-layer strategy ultimately achieves the same target rate-distortion point 
as sequential transmission, its progressive performance is superior. 

Although in this article we focus on the MSE distortion measure, image quality is subjective and 
depends on the application. Neither a thousand words nor a single scalar metric can accurately describe 
the “quality” of an image. 


III. Subband Coding 

A progressive quantization and compression scheme, such as the technique illustrated in the example 
of Section II, can be used in combination with a decorrelating transform. In particular, we use a subband 
coding stage prior to quantization and compression. Unlike block-transform-based methods, such as the 
discrete cosine transform used in the Joint Photographic Experts Group (JPEG) algorithm, subband 
coding does not suffer from blockiness artifacts even when used in progressive transmission. For a general 
description of image compression via subband coding, see [9]. 

A two-band subband decomposition uses high-pass and low-pass digital filters to decompose a data 
sequence into high and low subbands, each containing half as many points as the original sequence. This is 
done independently on horizontal and vertical lines of the image. Because image signal energy is usually 
concentrated in the lower frequencies, the lowest subband may be repeatedly decomposed. Figure 2 
illustrates a two-stage decomposition of the test image. 

The filters used for the tests in this article 1 are eighth-order quadrature mirror filters from [1, p. 267]. 
Because quadrature mirror filters are orthogonal, the MSE introduced by quantization in the transform 
domain is equal to the MSE of the reconstructed image. We can use this fact to improve the progressive 
rate-distortion performance by carefully selecting the order in which blocks of information from each 
subband are transmitted. At each step, we transmit the set of compressed bits giving the largest reduction 
in MSE per bit. Nonorthogonal filters have the potential to offer improved performance for the same 
complexity [7]; however, when such filters are used, the analysis in the subsequent sections no longer 
applies. 

The filters are implemented using circular convolution, i.e., each data block is periodically extended 
before filtering [6]. At the edges of data blocks, this often produces high-frequency components that are 
processed separately from the rest of the subband data because they are not as easily compressed. Since 
these components form a small fraction of the image, the penalty for inefficient compression is small. An 
elegant alternative to circular convolution is to use interpolation at the edges [10]. Since the improvement 
in rate-distortion performance is small, this alternative is not investigated in this article, though it could 
be implemented without significant change in the progressive transmission strategy described here. 

Each subband generally has a probability density function (PDF) quite close to Laplacian. 2 Figure 3 
illustrates an empirical PDF for one of the subbands of the test image. In the following sections, we will 
analyze a Laplacian source to determine a strategy for the transmission order of the subband data. 


1 For this article, we are more concerned with the overall compression and transmission strategy than with the selection of 
a particular filter. 

2 The lowest subband may be less well behaved and can often be compressed more efficiently after taking differences, 

although this process makes progressive transmission more difficult. Some performance improvement may be possible by 
revising the coding strategy used for the lowest subband. 
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Fig. 2. Two-stage decomposition: (a) original image and (b) image decomposed into 
subbands labeled a-g (contrast enhanced). 



TRANSFORMED VALUE 


Fig. 3. Empirical PDF of the subband 
labeled "b" in Fig. 2. 


IV. A Block-Oriented Progressive Transmission Scheme 

A simple technique for progressive transmission using subband coding is to transmit the subbands in 
order from lowest to highest frequency [9]. Moderate quality images can be reconstructed even when only 
the lowest subband has been received, because signal energy in images is usually more concentrated in 
the lower frequencies. However, the ability to transmit each subband progressively provides an added 
dimension that can improve the progressive rate-distortion performance because we can switch between 
the subbands during transmission. In this section, we describe a progressive transmission method for an 
individual subband. 

The progressive transmission scheme we use is bit-wise arithmetic coding [5]. A source sample is 
quantized, and each quantizer interval is identified by a 6-bit codeword. Figure 4 illustrates the codeword 
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Fig. 4. Codeword assignment for a 4- bit quantizer. 


assignment for a 4-bit uniform scalar quantizer along with the PDF f(x) for a Laplacian source with 
variance 2/a 2 . We can see that the bit assignment scheme is progressive — deleting i bits from the end 
of each codeword produces the same effect as using a lower resolution b — i bit quantizer. The codeword 
assignment scheme also ensures that, for a Laplacian source, a zero is more likely in every bit position, 
with the exception of the “sign” bit. 

The codewords corresponding to the quantized coefficients of a subband are grouped together. The 
zth bit layer, which consists of the zth bit from each codeword in the group, is compressed using the 
block-adaptive binary arithmetic encoder described in [5]. Each layer is compressed independently, i.e., 
the arithmetic encoder uses an estimate of the unconditional probability of a zero at each layer. 

Slight modification of the quantizer can improve the progressive rate-distortion performance. Figure 5 
shows three different quantizer options (illustrated for 3-bit quantizers). We can evaluate the rate and 
MSE distortion of each quantizer using the results of Appendix A. In each case, the quantizer range 
[— W, W], step size A, and number of bits b are related by 

A = W2 1-6 

A continuous rate-distortion performance curve for a quantizer can be obtained by fixing b and varying 
A. Figure 6 shows several such curves for the quantizer illustrated in Fig. 5(c). Decreasing A increases 
rate and lowers distortion until a minimum distortion point is reached. After this point, the distortion 
increases and the quantizer becomes inefficient. Ordinarily this inefficient region is not shown. 

When we use a quantizer progressively, the range is fixed, but the number of bits b increases, reducing 
step size by half for each successive bit layer transmitted. The points in Fig. 6 illustrate progressive rate- 
distortion performance. We observe in the figure that as b becomes large, the distortion does not approach 
zero, but is dominated by the contribution from the overload regions. As this happens, the progressive 
rate-distortion performance flattens, as points are lying along the inefficient regions of the rate-distortion 
curves. This floor is easy to find and depends only on the range and the Laplacian parameter a. We find 
from Eq. (A-l) in Appendix A that 


, MSE 
lim — 

6— KX> ( 7 * 
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Fig. 5. Different quantizer options: (a) uniform symmetric with bin boundary at origin, (b) uniform with reconstruction 
point at origin, and (c) almost uniform with enlarged center bin. 



when the reconstruction point of each interval is the interval midpoint. 3 In a real application, the range 
of possible values is limited by instrument dynamic range and the filters chosen, so the MSE floor can be 
avoided by making the quantizer range sufficiently large. 

Figure 7 compares the performance of each quantizer option of Fig. 5. The rate-distortion curves 
shown correspond to 8-bit quantizers. The progressive performance points are obtained by transmitting 
bit-layers with a target rate of 6 bits/sample. In Appendix A, we derive rate-distortion functions used to 
compare the quantizers. 

Quantizer (a), a symmetric uniform quantizer, has poor performance at low bit rates, because the 
sign bit is incompressible; hence, rates below 1 bit/sample are not achievable. Quantizer (b) is obtained 
by shifting quantizer (a) by A/2, so that a reconstruction point is present at the origin, as done in [5]. 
When A is large, low rates are achievable, and the performance obtained by varying A is close to the 
rate-distortion limit. However, when this quantizer is used progressively and the step size is small, the 
entropy of the sign bit approaches 1. Thus, the progressive performance is generally poor at low bit rates. 


3 This quantity is reduced by half when the centroid of each quantizer interval is used (see Eq. (A-2) in Appendix A). 
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Fig. 7. Rate-distortion performance of different 8-bit 
quantizers. The continuous curves are obtained by 
varying the quantizer step size A while keeping the 
number of bits fixed. The progressive performance 
points are obtained by calculating the rate-distortion 
performance after the transmission of each bit layer. 


To overcome these problems, we use quantizer (c) with range [— W, W], which effectively combines the 
two center regions of quantizer (a), in combination with a different transmission order. We do not initially 
transmit the sign bit layer, but rather we begin by transmitting the next layer (after arithmetic coding). 
Then, for each codeword where a “1” appeared (i.e., each quantized value in the range [-W, -W/ 2] U 
[W/2, IF]), the sign bit is transmitted. Then the next layer is transmitted, followed by the sign bits for 
each codeword representing a quantized value in the range [— W/2, —IF/ 4] U [IF/4, W/2\. This continues 
at each layer. Finally, if the quantized value lies in the range [— W2 1 ~ b , W2 1 ~ b ] (the two centermost 
quantization points), the sign bit is never transmitted. In progressive reconstruction, for any quantized 
value for which no sign bit was received, the origin is used as the reconstruction point. The advantage of 
this transmission order is that we delay the transmission of the sign bits, which are not easily compressible. 
Figure 7 shows that this technique offers improved progressive performance compared to quantizers (a) 
and (b) and is, in fact, close to the rate-distortion limit for a Laplacian source with MSE distortion, which 
is computed in Appendix B. 


V. Ordering the Compressed Subband Data 

In the previous section, we described a progressive transmission strategy that can be used for each 
subband. This strategy allows us to interrupt transmission of a subband to transmit data from another 
subband. In this section, we describe a method for choosing the order of transmission of subband data 
designed to improve the progressive rate-distortion performance. 

At any point in transmission, we wish to transmit the next (compressed) layer of bits from the subband 
that gives the largest reduction in MSE per transmitted bit. Thus, we select the bit layer from the 
subband with a rate-distortion curve whose slope is minimum. Using the analytical expressions for rate 
and distortion of the quantizer (see the Appendices) is rather intractable, so we use the approximation 

MSE ta a 2 e~ pR = \e~ pR 
a z 

where R is the rate in bits, and (3 = 1.3. This approximation is shown in Fig. 6. Using this approximation, 
omitting some algebraic details, we have the following subband selection strategy: Transmit the next bit 
layer from the subband with parameter a and rate R that minimizes In a + (3R/ 2. 
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To apply this strategy, we need to estimate the Laplacian parameter a for each subband or, equivalently, 
to estimate the mean absolute value 1/a. For each subband, we update the estimate at each bit layer 
transmitted. The estimate is created by keeping track of the frequency of sign bits transmitted at each 
bit layer, which is a quantity already computed in the compression stage. 



Fig. 8. Frequencies used to estimate the Laplacian 
parameter a at the third bit layer. 


Let fi denote the frequency of sign bits transmitted at the zth layer, and fc^ denote the frequency 
of the center region at the zth layer. This is illustrated in Fig. 8. For a Laplacian distribution with 
parameter a, at the zth bit layer, we expect 




-ocW/2 


3 = 1 


«W2-’ _ e -aW2'~^ 1C j<i 


/i° =1 — e 
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which gives parameter estimates of 
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(i) 


where, in Eq. (1), we choose the estimator closest to a^\ For our overall estimate of the mean absolute 
value at the zth bit layer, we take the weighted average of the above estimators: 


-1 J?(*) ^ £ 

i. _ h y ^ It 

“ a { c ] aj 


( 2 ) 


If any of the fi ’ s are zero, we can assign any finite value to the corresponding dj, since it will have no 
weight in Eq. (2). 


95 


Since the quantizer range may be quite large compared to the variance of a subband, it is common 
for the first few bit layers to consist entirely of zeros. For each subband, we first transmit an integer 
identifying the number of leading all-zero bit layers. This ensures that a is always well defined. 

VI. Results 

Curve (a) in Fig. 9 illustrates the progressive rate-distortion performance on the test image obtained 
using the subband selection strategy described in Section V combined with the parameter estimates from 
Eq. (2). In the figure, we compare this performance to the more traditional progressive transmission 
method (transmitting the subbands sequentially from lowest to highest frequency), shown as curve (b), 
and the optimum transmission order (that obtained by testing every possible transmission order of the 
subband data), curve (c). We also show the performance of the graphic-in-line-format (GIF) compression 
technique, a progressive transmission method that has recently become popular in Web browsing software. 

In Fig. 10, we show four intermediate reconstructed images obtained in progressive transmission. The 
rate-distortion points corresponding to these four images are shown in Fig. 9. 



Fig. 9. Progressive rate-distortion performance on the test image using 
12-bit quantizers on each subband. Curve (a) shows the performance 
obtained using the strategy described in Section V. Curve (b) is the 
performance obtained when we transmit the subbands in order from the 
lowest to the highest frequency. Curve (c) is the performance obtained if 
we choose the optimum subband bit layer at each step. The curve labled 
"GIF" corresponds to the performance obtained using the GIF progressive 
transmission technique. The images corresponding to the four labled 
reconstruction points are shown in Fig. 10. 


VII. Conclusion 

We have shown that careful selection of the quantization and compression strategy can significantly 
improve the progressive rate-distortion performance of a data compression system. The use of transforms 
such as subband coding can lead to further improvements, because they provide additional freedom in 
choosing the order in which the compressed data are transmitted. 

When an effective progressive transmission scheme is used, the progressive rate-distortion performance 
may be good at every point in transmission. For image compression, this means that even if the rate 
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Fig. 10. Reconstructed images at various stages in the transmission. The savings compared to transmission of 
the original 8-bit image without compression are (a) 18.2, (b) 16.5, (c) 15.5, and (d) 13.9 dB, respectively. 


constraints change during a mission, redesign of the compression strategy is not necessary. When the 
reverse channel can be used for browsing, the benefits of progressive transmission become even greater — 
we are not required to decide in advance the amount of resources devoted to a particular image. 

It is interesting to consider the savings in dB that can be offered by lossy data compression. Suppose 
that 2:1 lossless compression is possible for the original image shown in Fig. 2. In this case, the lower- 
quality image of Fig. 10(d) could be obtained at a savings of approximately 10.9 dB. Alternatively, for 
the price of a single losslessly compressed image, we could obtain 12 images of quality comparable to that 
of Fig. 10(d). 
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Appendix A 

Performance of a Scalar Quantizer on a Laplacian Source 


In this appendix, we analyze the rate-distortion performance of a Laplacian source with PDF f(x) 
using a scalar quantizer that is nearly uniform. Consider the symmetric 6-bit quantizer (with 2 b — 1 
intervals) shown in Fig. A-l. The reconstruction point of each interval may be the midpoint or the 
centroid of the interval. 



Fig. A-1. A symmetric scalar quantizer. 


The normalized MSE obtained when this quantizer is used is 


MSE 


a 


2 


2 fc_l -3 

MSE c + 2MSE ov + 2 ]T MSEi 

i = 0 


where MSE C is the MSE contribution of the center region, 


MSE C = ~ [2 - (y 2 + 2 y + 2)] 


and y = aA/2. MSEi is the MSE contribution of the region [(w/2 + i) A, (rn/2 + i + 1)A), i = 
0, 1, • ■ • ,2 6_1 — 3, 


MSEi = 


o-y-si 


2a 2 


when the centroid is used, and 

-y-si 

MSEi = [(8 + s 2 )(l - e- s ) - 4s(l + e" 5 )] 


when the midpoint is used, where s = aA. MSE ov is the MSE of each overload region ±\(w/2 
+ 2 6 - 1 -2)A,oo), 
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MSE ov 


e s{2-M) e -y 

2ofi 


when the centroid is used, and 


MSE ov = 



-4s + 8\ 
8a 2 ) 


e s(2-2”-') e - y 


when the midpoint is used. 

Combining, we find that when midpoints are used, 


MSE 


1 — e~ y 


iy(y + 2) + ^s(4-s) + 


s 1 


s(2-2 l 




and when centroids are used, 


MSE 


1 

S 2 e s(' 1 _ e ,( 2 - 2b - i )y 

1 - -e~ y 
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{y ) + (e s — l) 2 


(A-l) 


(A-2) 


Next we compute the rate. The center region has probability P c — 1 — e~ y \ each overload region 
has probability P ov = (l/2)e~ y e s< - 2 ~ 2b ~^; and the region \{w/2 + z)A, (w/2 + i + 1)A) has probability 
{l/2)e~ y ~ i3 {\ - e -s ). The entropy is 


H = ( e~ y - 1) ln(l - e- y ) + e~ y {y + In 2) + e~ y (l - e^ 2 " 2 ' ')) 



— ln(l — e s ) 


nats 


Using the bit-wise transmission strategy described in Section IV, the sign bit has entropy of 1 bit but 
is not transmitted when the sample value lies in the center region; this contributes e~ y bits to the rate. 
If j denotes the bit layer index, with j — 1 corresponding to the least significant layer and j = b the sign 
bit layer, then the probability that the jth bit is a 1 is 


so the rate of bit-wise transmission is 


6-2 

R = e~ y + Y,'H2 

j = o 


-s2 j 


o a ~V 


+ e 


-s 2° 


1 + e _s23 


bits 


where Ti 2 ( x ) = 


— log 2 x — (1 — x) log 2 (l — x ) is the binary entropy function. 
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These expressions are used to compare the different quantizers in Section IV. We can also observe the 
losses resulting from the bit-wise transmission strategy and using the midpoint rather than the centroid of 
each quantizer interval. In fact, these losses are negligible, as can be seen by comparing the rate-distortion 
curves in Fig. 7. 


Appendix B 

Rate-Distortion Limit For a Laplacian Source With MSE Distortion 


In this appendix, we derive the rate-distortion limit of a Laplacian source using a scalar quantizer with 
MSE distortion metric. A Laplacian source with parameter a has PDF f(x) = (a/2)e~ a ^ and variance 
a 2 = 2/a 2 . For a discussion of the mean absolute error distortion metric for this source, see [3]. 

The goal is to find the optimum (possibly infinite) scalar quantizer. Let the quantizer regions be 
[L,L+i), *=••', -1,0, 1,2, • • •. The MSE distortion of the quantizer is MSE = Y^Z-oo MSE it where 
MSEi is the contribution to the MSE from the region [t t , fj + i) (assume 0 <U < U+i] a similar analysis 
applies when this is not the case): 


MSEi 


fU+i 
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(x - Ci) 2 f(x)dx 
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Here c* is the centroid of the region (this choice minimizes MSE): 
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and Pi = (1/2) (e ati — e ati+1 ) is the probability of the region. The total entropy associated with the 
quantizer is If = — na -ts- 

Consider what happens when we vary the threshold f»: 
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So by the chain rule, the slope of the rate-distortion curve (when varying £» ) is dMSE/dH = 
( dMSE/dti)/dH/dti . If the quantizer is locally uniform with step size A, i.e., fj_i =U — A, t i+ 1 — U+A, 
then (omitting algebraic manipulations), 


dMSE (2 - s)e s -2-s 
dH ~ a 2 (e s - 1) 


where s = aA. By the Kuhn-Tucker conditions, since this expression depends only on the normalized 
step size s, the optimal scalar quantizer is uniform (and infinite) on both sides of the origin and has the 
same step size on both sides. 

Next, we need to determine the behavior of the quantizer near the origin. Suppose there is an arbitrary 
(possibly different) offset on each side of the origin. By symmetry, setting these two offsets equal will give 
a solution to the Kuhn-Tucker conditions. So we want to know the width of the center region relative to 
the normalized step size s. 

Let the width of the center region be wA. The center region has probability mass 1 — e~ y , where 
y = sw/ 2, and makes MSE contribution 

/•tuA/2 i 

MSE C = 2 J x 2 f(x)dx = -[2-(y 2 + 2y + 2) e~ y ] 

The right-hand side regions are [(w/2 + i) A, (w/2 + i + I)A), i = 0, 1, 2, • • •. Using i* = (to/2 + i) A and 
ti+i = (w/2 + i + 1)A in Eq. (B-l) yields 


MSE, 


o-si-y 

2o? 


1 — e~ s + 


1 — e s 


Pi= 2 e ~ Si ~ y ( 1 ~ e ~ S ) 


so the rate contribution is Ri = -Pi In Pi = (l/2)e 31 y (l-e s ) [y + si + ln2 - ln(l - e~ fl )]. The total 
MSE is 


OO 

MSE = MSE C + 2^2 MSEi 
i = o 


M 2+e -’ 


L(1 - e s )(l - e~ 5 ) 



(Alternatively, this can be derived from Eq. (A-2) in Appendix A taking the limit as 6 — > oo.) The 
entropy is 
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oo , 

-Pc In P c + 2 Y,(~ p i ln p i) =e ~ V \- { 

i=0 ^ 


+ In 


f2.(e y -l) 


- ln(l - e- y ) 


Given these expressions for MSE and entropy in terms of s and y, we apply the Kuhn-Tucker conditions 
once more. When s is fixed and y varies, the resulting rate-distortion curve has slope 


102 



s 2 e s 


(e s - l) 2 


e s - 1 


+ In 


2(e y — 1) 


1 - e _s 


When y is fixed and s varies, the resulting rate-distortion curve has slope 



By the Kuhn-Tucker conditions, these two quantities must be equal. We can use this fact to solve for y 
in terms of s numerically. We find that w varies between 1 (high rate limit) and 2 (low rate limit). The 
resulting rate-distortion limit is shown in Fig. 7. 
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Ka-Band Monopulse Antenna-Pointing Systems 
Analysis and Simulation 

V. Y. Lo 

Communications Systems and Research Section 


NASA ’s Deep Space Network (DSN) has been using both 70-m and 34-m reflec- 
tor antennas to communicate with spacecraft at S-band (2.3 GHz) and X-band 
(8.45 GHz). To improve the quality of telecommunication and to meet future 
mission requirements, JPL has been developing 34-m Ka-band (32-GHz) beam- 
waveguide antennas. Presently, antenna pointing operates in either the open-loop 
mode with blind pointing using navigation predicts or the closed-loop mode with 
conical scan (conscan). Pointing accuracy under normal conscan operating condi- 
tions is in the neighborhood of 5 mdeg. This is acceptable at S- and X-bands, but 
not enough at Ka-band. Due to the narrow beamwidth at Ka-band, it is important 
to improve pointing accuracy significantly (~2 mdeg). Monopulse antenna tracking 
is one scheme being developed to meet the stringent pointing-accuracy requirement 
at Ka-band. Other advantages of monopulse tracking include low sensitivity to 
signal amplitude fluctuations as well as single-pulse processing for acquisition and 
tracking. This article presents system modeling, signal processing, simulation, and 
implementation of Ka-band monopulse tracking feed for antennas in NASA/DSN 
ground stations. 


I. Introduction 

The design of the DSN monopulse pointing system consists of the reflector antenna, multimode corru- 
gated horn feed, waveguide coupler, monopulse signal processor, and other associated RF electronics. A 
general block diagram is provided in Fig. 1. Starting at the main reflector, a tapered beam is formed. The 
HE11 mode in the corrugated horn is excited to radiate the sum pattern while the TE21 mode waveguide 
coupler generates the difference pattern [1]. With the assumption of perfect Ka-band-to-IF conversion, 
signal processing starts in the IF domain. A phase-locked loop recovers the carrier phase. This is used 
as a phase reference to coherently demodulate the elevation and azimuth difference channels. The sum 
and difference baseband signals are low-pass filtered and then followed by monopulse signal processing 
from which elevation and azimuth pointing errors are estimated. The error signals are used to drive the 
antenna servo controller for pointing corrections. 

To predict the overall system performance, an antenna system model is needed. Mathematical models 
of the open-loop S-curve and pointing variance are presented in Section II. These models are validated 
by comparison with the difference pattern generated by a single-aperture multimode antenna. In Sec- 
tion III, a simplified servo loop model on the antenna controller and driver system is described. Higher- 
order complex multistate models are reserved for future upgrade. An end-to-end DSN block simulation 
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Fig. 1. Antenna-pointing system with monopulse signal processing. 


program is presented in Section IV. The system includes the monopulse antenna block for a single-aperture 
multimode antenna, the digital receiver block for the Block V receiver (BVR), the decoder block for the 
maximum-likelihood convolutional decoder (MCD), the low-pass filter (LPF) and the monopulse processor 
blocks for pointing-error estimation, and the servo block for the antenna driver controller. Finally, in 
Section V, open- and closed-loop system performance in the presence of noise, cross-channel interference, 
amplitude and phase imbalance, and wind loading are investigated using this block simulation program. 


il. System Model of Monopulse Antenna Pointing 

Based on the physics of a corrugated horn and mode coupler, the classical four-horn monopulse antenna 
is shown to produce the same open-loop S-curve as the monopulse single-aperture multimode antenna. 1 
In the presence of random noise, it can also be shown that the carrier-to-noise power spectral density ratio 
(CNR) is identical between the four-horn model and the single-aperture model with the same aperture 
size. Thus, the four-horn model produces the same system performance as a single-aperture multimode 
antenna under identical input signal conditions. The simple four-horn model can now be used for system 
simulation and analysis. Without loss of generality, we shall examine the azimuth pointing errors only. 

From a complex signal analysis on a four-horn model [2, 3], 2 the mean error voltage and its variance 
are shown to be 


1 S-curves from various models were compared under deterministic conditions in T. K. Wu, “Monopulse Antenna Study,” 
JPL Interoffice Memorandum 3365-94-MP-001 (internal document), Jet Propulsion Laboratory, Pasadena, California, 
September 1994, and in V. Lo, “Single Aperture Multi-Mode Monopulse Antenna Pointing With Corrugated Horn,” 
JPL Interoffice Memorandum 3393-94- VYL-015 (internal document), Jet Propulsion Laboratory, Pasadena, California, 
November 1994. 

2 The detailed mathematical derivation is shown in V. Y. Lo, “Monopulse Tracking Feed Link Study,” JPL Interoffice 
Memorandum 339-92-134 (internal document), Jet Propulsion Laboratory, Pasadena, California, December 1992. 
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pointing error, the pointing error and its variance simplify to 
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The 34-m Ka-band beam-waveguide (BWG) antenna parameters can now be substituted into Eqs. (1) 
and (2). For d/X = 1000, the pointing-error voltage plot from Eq. (la) is compared with S-curves 
generated from other physics and SPW 3 simulation models in Fig. 2. It shows that, at medium pointing 
errors, the mathematical model from Eq. (la) matches the physics model of a single-aperture multimode 
antenna as well as the SPW-simulated S-curve of a four-horn monopulse antenna. Since a monopulse 
processor essentially performs integration on the input signal CNR to optimized pointing-error estimates, 
Eq. (2b) can be rewritten in terms of integration time T, ag = (18xl0~ 3 ) \\/ \J{CNR ■ T) . The linearized 
standard deviation is plotted against the integration time for CNR values ranging from 10 to 30 dB-Hz 


1 Signal Processing WorkSystem (SPW) is the registered trademark of Alta Group, formerly Comdisco Systems, Inc. 
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(see Fig. 3). Based on the DSS-13 KABLE 4 link prediction, the carrier-to-noise power spectral-density 
ratio is expected to be greater than 10 dB-Hz. For a CNR value of 10 dB-Hz and a pointing error of 
2 mdeg, integration time is 8 s. Given the possible range of CNR, 0.08 to 8 s represents the corresponding 
spread in integration time. For stable loop operation, the observation time of the estimator has to be 
short compared to the update time of the servo loop. 


i i I i i i i i i i i i | i i i i 

■ ANTENNA POINTING SIMULATION (SPW) 

■ SMALL POINTING-ERROR MODEL (FOUR-HORN) 

1- -£r MEDIUM POINTING-ERROR MODEL (FOUR-HORN) 

■ SINGLE-APERTURE MULTIMODE MODEL 
(CORRUGATED HORN AND REFLECTOR) 


i | i i i i i i i i i 



POINTING ERROR, deg 


Fig. 2. Comparison of simulation, small/medium pointing-error 
models, and the single-aperture multimode model. 



III. Simplified Servo Loop Model of the DSN Ground Antenna 

A model of the DSN BWG antenna servo control system [4] based on DSS 13 was developed to study 
pointing dynamics, channel cross-coupled dynamics, and wind disturbances. The state-space model of the 
antenna structure was obtained from its finite-element model. State reduction was applied separately on 
the antenna structure, elevation and azimuth drives, and rate-loop model. To a first-order approximation, 


4 The Ka-band link experiment (KABLE) was presented in OSC Advanced Systems Review: TDA Systems Development 
(internal document), Jet Propulsion Laboratory, Pasadena, California, June 1992. 
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a proportional controller is used as a pointing-error estimator based on the observed differential error 
voltages. Such information is passed on to the elevation and azimuth drives to turn the antenna structure. 
This second-order servo loop model is shown in Fig. 4. The frequency response has a 3-dB roll-off at 
0.2 Hz. 


IV. End-to-End DSN Ground-Link Simulation 

APSIM_MON is a simulation software [5] of the DSN BWG antenna supporting the Ka-band monopulse 
development effort. It is constructed on the platform of the SPW. The hierarchy simulation model for 
monopulse antenna pointing (ANT.MON) has been integrated with the advanced digital receiver model 
(DIG-RCVR) and the decoder model (CONJDEC) to form an end-to-end DSN telemetry system (see 
Fig. 5). 5 The proper encoded signal is available through the simulation signal generator (CC-SIG) model, 
which generates an encoded pseudonoise (PN) data sequence with different initial phase and frequency 
conditions. Parameters of the SNR, loop acquisition, and running averages of the bit error rate are 
monitored by the simulation monitor block (PERF_MONITOR). 

The simulation model has been applied to analyze various design parameters. The following cases have 
specifically been investigated: 


(1) Open- loop monopulse antenna pointing with and without random noise 

(2) Closed-loop monopulse antenna pointing with and without random noise 

(3) Effects of amplitude and phase imbalance between the sum and difference channels 

(4) Effects of pointing errors on end-to-end telemetry link performance 


5 The SPW Costas-loop digital receiver model developed by J. Gevargiz for the BVR has been expanded to include residue 
carrier with square- and sine-wave subcarriers; see J. Gevargiz, “Acquisition, Tracking and Bit-Error-Rate Analysis of the 
Block V Receiver’s BPSK Tracking Loop Including the Convolutional Coder and Decoder,” JPL Interoffice Memorandum 
3396-93-05 (internal document), Jet Propulsion Laboratory, Pasadena, California, March 1993, and V. Lo, “Antenna 
Pointing Simulation With Monopulse Processing (APSIMJvlON),” JPL Interoffice Memorandum 3393-95- VYL-02 (inter- 
nal document), Jet Propulsion Laboratory, Pasadena, California, April 1995. 
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Fig. 5. SPW simulation of the monopulse antenna-pointing system. 


V. Summary of Simulation Results 

For case (1), with variable pointing errors, the corresponding error voltages in the forward loop trace 
out the S-curve shown in Fig. 2. When Gaussian noise is added into the loop, a noisy pointing-error 
voltage is observed. The standard deviation of the linearized mathematical prediction based on the 
corrupt pointing-error voltage, Eq. (2), is shown in Fig. 3. Knowing the input SNR, these curves provide 
the integration time needed to meet a specific pointing-error variance. 

For case (2), a simple second-order servo loop is used to model the antenna controller. The frequency 
response has a 3-dB roll-off at 0.2 Hz (see Fig. 4). Stability of the servo loop is investigated as a 
function of the low-pass filter bandwidths in the forward path. Various types of low-pass filters, including 
Butterworth, Bessel, and Chebyshev, with up to twenty orders are used in this evaluation. Deterministic 
step responses of the loop are simulated (see Fig. 6). The results show that the low-pass filter bandwidth 
has to be about two orders of magnitude higher than the servo-loop bandwidth for stable loop operations. 
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Fig. 7. Standard deviation of closed-loop pointing error of the azimuth 
channel induced by step perturbation in single and dual channels. 


The same conclusion can also be reached from root locus and Nyquist stability analyses. In the presence 
of uncorrelated Gaussian noise, the simulated standard deviations of closed-loop step pointing error for 
single-channel and cross-channel interference are shown in Fig. 7. Coupling between the azimuth and 
elevation channels will increase the pointing error beyond this baseline. 

For case (3), the amplitude imbalance between the sum and difference channels scales the voltage of 
the S-curve while the relative phase imbalance in electrical length manifests as a bias of the pointing error. 
The source of such imbalance may be due to amplifiers, hybrids, couplers, filters, etc. Both open- and 
closed-loop pointing-error variances are also modified by a multiplication factor of the squared amplitude 
imbalance. These results are summarized in Fig. 8, where a relative amplitude of 1 and a relative phase 
of 0 deg represent the baseline with no imbalance. Offsets corresponding to 0.1- and 1-dB pointing losses 
are marked by dotted lines in the figure for reference. Specifications on amplitude and phase mismatch 
between the difference and sum channels can now be established to meet pointing-loss tolerance. 
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Fig. 9. System loss due to step pointing error. 

Finally, the impact of pointing errors on the link quality driven by wind loading or gravitational effects 
can now be analyzed by observing the relative increase in bit errors after the onset of such disturbances. 
The loss in signal power relative to the unperturbed case can be obtained from the increase in bit-error 
probability. For example, a 1-kbps downlink with a (7,1/2) convolutional code at an 80-deg modulation 
index in a monopulse closed-loop pointing system with system loss from step pointing error induced by 
wind load is shown in Fig. 9. Unlike the antenna pointing loss from miscalibration, this peak pointing 
loss represents a dynamic process depending on the rate of change of the disturbance relative to the servo 
loop tracking rate. Similar to radio loss in carrier tracking, the results show that the telemetry link is 
more susceptible to small pointing errors at high SNR. 


Ill 



VI. Conclusion 


A system model for monopulse antenna pointing has been presented. Under the framework of SPW 
simulation, a monopulse antenna block is constructed from this system model. Both open- and closed-loop 
responses in the presence of noise are investigated through simulation. Pointing performance is evaluated 
from the perspective of SNR degradation (system loss) on an end-to-end link due to wind loading. These 
simulation results provide the basis for design, implementation, and performance evaluation of the DSN 
Ka-band monopulse pointing system. 


Acknowledgments 

I would like to acknowledge the technical support by M. K. Sue and T. K. Wu. 


References 


[1] P. J. Clarricoats and A. D. Olver, Corrugated Homs for Microwave Antennas, 
London: Peter Peregrinus Ltd., 1984. 

[2] T. Inoue and T. Kaitsuka, “K-Band Tracking System for Domestic Satellite 
Communication System,” IEEE Trans, on Aerospace and Electronic Systems, 
vol. AES-17, no. 4, pp. 561-570, July 1981. 

[3] S. Sharenson, “Angle Estimation Accuracy With a Monopulse Radar in the 
Search Mode,” IRE Trans, on Aerospace and Navigation Electronics, vol. ANE9, 
no. 3, pp. 175-179, September 1962. 

[4] W. Gawronski and J. A. Mellstrom, “Modeling and Simulations of the DSS 13 
Antenna Control System,” The Telecommunications and Data Acquisition Prog- 
ress Report 4.2-106, April- June 1991 , Jet Propulsion Laboratory, Pasadena, Cal- 
ifornia, pp. 205-248, August 15, 1991. 

[5] V. Y. Lo and M. K. Sue, “Monopulse Signal Processing and Simulation for 
DSN Beam Waveguide Antenna,” Proc. ICSPAT, vol. 2, Boston, Massachusetts, 
pp. 1603-1607, October 1995. 


112 



TDA Progress Report 42-124 



a 


/ 


February 15, 1996 


Modeling and Analysis of the DSS-14 
Antenna Control System 

W. Gawronski and R. Bartos 
Communications Ground Systems Section 


An improvement of pointing precision of the DSS-14 antenna is planned for the 
near future. In order to analyze the improvement limits and to design new con- 
trollers, a precise model of the antenna and the servo is developed, including a finite 
element model of the antenna structure and detailed models of the hydraulic drives 
and electronic parts. The DSS-14 antenna control system has two modes of opera- 
tion: computer mode and precision mode. The principal goal of this investigation 
is to develop the model of the computer mode and to evaluate its performance. 
The DSS-14 antenna computer model consists of the antenna structure and drives 
in azimuth and elevation. For this model, the position servo loop is derived, and 
simulations of the closed-loop antenna dynamics are presented. The model is sig- 
nificantly different from that for the 34-m beam-waveguide antennas. 


I. Introduction 

The DSS-14 antenna control system model consists of the antenna structure, antenna drives in azimuth 
and elevation, and the position servo loop. Each drive, in turn, consists of gearboxes, hydraulic servo 
(active and passive valves, hydraulic lines, and hydraulic motors), and electronics boards (amplifiers 
and filters). The DSS-14 antenna control system model was developed by R. E. Hill [1,2]. In the present 
development, we obtain a more precise model that allows for accurate simulations of the antenna pointing 
errors and allows simulation of the intermediate variables, such as torques, currents, wheel rates, truss 
stresses, etc. We incorporate the finite element structural model with free rotation in azimuth and 
elevation, in a manner similar to the 34-m antenna models [3-5], that involves cross-coupling effects 
between azimuth and elevation, wind pressure on the dish, and pointing error model. The hydraulic part 
involves a recent development in modeling of the hydraulic components by R. Bartos [6-8]. 

The rate loop model consists of the elevation and azimuth drives and the antenna structure. Each drive 
consists of three major components: the electronics boards, hydraulic system, and gearbox. A model of 
each component is derived separately, then put together, forming the drive and rate loop models. Finally, 
the position loop is closed to obtain the position loop model. 


II. Drive Model 

A block diagram of the drive model is shown in Fig. 1, where N t is the ratio between motor rate 
and tachometer (pinion) rate; r, rad/s, is the rate input to the drive; i, A, is the hydraulic active valve 
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Fig. 1 . Block diagram of the antenna drive. 


solenoid current; T 0 and T, N-m or lb in., are the gearbox and on-axis torques, respectively; and d tac h and 
6 m , rad/s, are the pinion and motor rates, respectively. The state-space representation of the electronic 
board, hydraulic system, and gearbox are derived in the following sections. 

A. Electronic Board 

A schematic diagram for the electronic board is shown in Fig. 2. The inputs are the rate command 
r, rad/s, and the tachometer rate 6 tac h , rad/s. The output is the solenoid valve current i. The scaling 
factors, k r and k t , convert the inputs into the command voltage, v r , and tachometer voltage, v t . The 
subsystem, G t , is the tachometer circuit: it transforms the tachometer voltage, vt, into the voltage, v t0 . 
The subsystems with the transfer functions G r i and G r 2 are the rate amplifier circuits: they transform 
the command voltage, v r , and the tachometer voltage, v t0 , into the error voltage, v s . The subsystem with 
the transfer function, G s , is the valve driver amplifier circuit, with the error voltage, v 3 , as the input and 
the valve current, i, as the output. 



Fig. 2. Block diagram of the electronic board. 


The following transfer functions of each of the four components are derived in the Appendix. The 
transfer function Gto for azimuth is 


G to = 0.151 


( 1 ) 


and, for elevation, is 


G to = 0.127 


( 2 ) 


The transfer functions G rl (from v r to v s ) and G r2 (from v to to v 3 ) are 
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G r 1 = 6.20 Go ' 
G r2 = — 4.65 G 0 t 


where 


1 4- 0.400s 
1 4- 4.205s 


is the transfer function of a lag compensator. The transfer function G a is 


G s = 4.42 x 10~ 5 


( 3 ) 


( 4 ) 


( 5 ) 


The scaling factors, k r and k t , are k r = 1212.6 V/rad/s and kt = 2.5 V/rad/s. Thus, the command 
transfer function from the rate command r to the solenoid current i s is 

G r = k r G rl G s = 0.3323 G 0 (6) 

where G 0 is defined in Eq. (4). The tachometer transfer function from the tachometer rate Otach to the 
solenoid current i s is 


G t = k t G to G r2 G s (7) 

q _ ( —0.7750 x 10~ 4 G o for azimuth /g.. 

4 \ —0.6525 x 10~ 4 G o for elevation 

In order to check the correctness of the derivation, note that the ratio G r /G t should be equal to N 0 , where 
N 0 is the tachometer-to-axis ratio ( N 0 = 4287.5 for azimuth and N a = 5083.6 for elevation). Indeed, 
from Eqs. (6) and (7), one obtains 


G r ( —4287.7 = N 0 for azimuth 

G t \ —5092.7 = N 0 for elevation 


Finally, the state-space representations of the transfer functions G r and Gt (for azimuth and elevation) 
are easily obtained with the standard Matlab command in the form 


Xfo — AfoXb 4 “ B V 4 “ Bb2@tach 

> 

i = GfoX 4 “ T Db 2 6tach t 


( 10 ) 


The plot of the transfer function in azimuth (magnitude and phase) from r to i is shown in Fig. 3. The 
transfer function for elevation is identical. The plots of the transfer functions in azimuth and elevation 
from Otach to i are shown in Fig. 4. The magnitudes drop in the frequency range from 0.01 to 0.1 Hz due 
to implementation of the filter G c . 
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Fig. 3. Magnitude of the electronic board transfer function from the rate input to the 

board current. 



B. Hydraulic System 

The hydraulic servo system model was presented by R. E. Hill in [1] and [2]. Here we take a different 
approach, based on the recent investigations of hydraulic components by R. Bartos (see [6-8]). A block 
diagram of the DSS-14 hydraulic system is shown in Fig. 5. It consists of the hydraulic motor, shorting 
valve, hydraulic lines A and B, passive servo valves, and active servo valves. It has two inputs, servo 
valve current i and motor rate 6 m , and one output, motor torque T a . The equations for each component 
are derived separately based on the work of Bartos [6-8]. Basically, these models are nonlinear ones; 
however, we linearize them in order to model the antenna linear regime of operation. 

1. Active Servo Valve. This valve model has the input, i, A, and two outputs, q av — the flow rate 
out of port a, cm 1 * 3 /s, or in. 3 /s, and qb v , the flow rate out of port b, cm 3 /s, or in. 3 /s. From [6], one obtains 

Qav ‘^Co^oQav T U 0 ([ av — Ld 0 k a i (11a) 

Qbv “ Qav (lib) 


where ( 0 = 0.8 is the damping ratio, ui 0 = 345.6 rad/s is the valve natural frequency, and k a — 59,200 
— 97,300 cm 3 /s/A (23,300-38300 in. 3 /s/A) is the valve gain. The lower value is the Bartos estimate, 
while the upper value is the Hill estimate [2]. The values of the parameters are listed in Table 1. 
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Fig. 5. Block diagram of the hydraulic drive. 


Table 1. Parameters of the active servo valve (line A and line B). 


Drive 

Co 

O>o, 

rad/s 

fta.5 

cm 3 /s/A 

in. 3 /s/A 

Azimuth 

0.6 to 0.8 

345.56 

59,200 to 97,300 

23,300 to 38,300 

Elevation 

0.6 to 0.8 

345.56 

59,200 to 97,300 

23,300 to 38,300 


Introducing the new variable q v = q av — qb v , one obtains from Eq. (11) 

q v + 2( 0 u 0 q v -I ~ = 2w 0 fc 0 t (12) 

The differential variable q v and the other differential variables introduced allow one to further simplify 
the analysis without loss of accuracy and to get rid of the “parasitic” variables, such as tank pressure, 
supply pressure, and case pressure. 

2. Shorting Valve. The pressures p a and pb , kPa (lb/in. 2 ), are the inputs to the shorting valve, and 
the flows q as and qbs, cm 3 /s (in. 3 /s), are its output (see Fig. 5). The linearized relationship between the 
inputs and outputs is as follows: 


Qas — ks(Pa Pb) 1 


Qbs — Qas 


(13) 


where k s is the valve gain, k s = 0.0007 — 0.007 cm 3 /s/kPa (0.0003-0.003 in. 3 /s/psi), both in azimuth and 
elevation. 
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Introducing new differential variables p = p a ~Pb and q s = q as — qbs, one obtains Eq. (13) in the form 


q s = 2k s p 


(14) 


3. Passive Servo Valve. This valve has four inputs: pressures p a and pb, supply pressure p s , and 
tank pressure p t . The last two are supplementary constant inputs that can be removed from the analysis. 
The valve has two outputs: flows q ap and qb p . Its linearized input-output relationship is as follows: 


Qap = k p i(p a - p s ) + k p2 {p a - Pt) (15a) 

qbp = k p i(p b - p s ) + k p2 (pb - Pt) (15b) 


where the gains are k p \ = 0.0055 cm 3 /kPa (0.00233 in. 3 /psi) and k p2 = k p i, the supply pressure is 
17,240 kPa (2500 psi), and the tank pressure is 345 kPa (50 psi). These values are identical for azimuth 
and elevation. 

Introducing q p = q ap — qbp , and recalling that p = p a ~Pb , one obtains Eqs. (15a) and (15b) as follows: 


q P = (k p i + k p2 )p — 2 k p p 


(16) 


where, for simplicity of notation, we denote k p — k p i — k p 2 ■ 

4. Hydraulic Motor. The motor is described in [8]. From Fig. 5, it follows that the motor has four 
inputs and four outputs. The inputs are pressures p a and pb, case pressure p c , and motor rate 6 m , rad/s. 
The outputs are flows q a and qb, leakage to the case q c , and motor torque T 0 , N-m (or lb in.). Following 
[8], one obtains the flow q a from Eq. (59) of [8]: 


q a = q a i + q a 2 + q a 3 


(17) 


But, from Eq. (40) of [8], 


Qal — DQ m 


(18) 


where D = 6.3 cm 3 /rad (0.3836 in. 3 /rad) is the motor “displacement.” From Eq. (52) of [8], one obtains 


?a2 = ~~~~{Pa ~ Pc) (19) 

A 4 

where k a2 = from 6.35 x 10 -5 to 18.3 x 10 -5 cm 3 (from 2.5 x 10 -5 to 7.2 x 10 -5 in. 3 ) is the leakage 
constant (assumed to be 10 -4 cm 3 , or 4 x 10 -5 in. 3 ); p — from 2.8 x 10 -4 to 2.8 x 10~ 3 kPa s (from 
4 x 10 -5 to 4 x 10 -6 lb s/in. 2 ) is the absolute fluid viscosity (assumed to be 10 -4 cm 3 , or 4 x 10 -5 in. 3 ); 
and p c = 221 kPa (32 psi). From Eq. (58) of [8], one obtains the leakage from port A to port B, q a 3: 


ka3 ( , 

qa 3 = \Pa Pb) 

A* 


(20) 
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where k a 3 is the constant of proportionality determined through experiments. It is assumed to be equal 
to fc a 2 ) k/i '3 = k a 2 * 


Combining Eqs. (17) through (20), one obtains 


i k a 2 T kaZ k a 3 k a 2 

q a = D6 m + p a p b p c 

p p p 


( 21 ) 


From [8], Eq. (60), one obtains the flow rate <?{,: 


96 = 9bl + 962 “ 9 o3 


( 22 ) 


It follows from [8], Eq. (41), that 


96i — 9ai — D9 n 


(23) 


and from [8], Eq. (53), that 


962 = —(Pb-Pc) 
t 1 


(24) 


Combining Eqs. (22), (23), (24), and (20), one obtains 


A k a 3 k a 3 'r kb2 h,2 

q b = -DO p a + pb p c 

p p p 


(25) 


The motor torque T a is obtained from Eq. (28) of [8] by neglecting Coulomb friction and inertia torques 
(the latter are included in the gearbox model): 


T 0 = T P + T f (26) 

where T p is the torque generated by the motor and Tj is the viscous friction torque. The linearized 
Eq. (10) of [8] gives the torque generated by the motor: 

Tp — D(Pa ~ Pb) ( 27 ) 

and from Eq. (24) of [8], one obtains the viscous friction torque: 

T f = —k v Dp9 (28) 


where k v — 0.0438 is a dimensionless viscous friction coefficient. Combining Eqs. (26) through (28), one 
obtains 


T 0 = Dp a - Dp b - k v Dp9 


(29) 
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Define q = q a — qb\ then from Eqs. (21) and (25), one obtains 


o no i *^<*2 

q = 2 D9 m + p 

M 

T 0 = Dp - k v Dp9 m 


( 30 ) 


The motor parameters are given in Table 2. 


Table 2. Parameters of the hydraulic motor (line A and line B). 


Drive 

fi, kPa s 

H, lb s/in. 2 

D, 

cm 3 /rad 

D, 

in. 3 /rad 

k-y 

fc Q 2 , cm 3 

k a 2 , in. 3 

Pc, 

kPa 

Pci 

psi 

Azimuth 

2.8 X 10~ 4 

0.4 x 10~ 6 

25.2 

1.52 

0.0438 

4.1 x 10“ 4 

2.5 x 10- 5 

220 

32 


to 

to 




to 

to 




28 x 10" 4 

4 x 10" 6 




11.8 x 10“ 4 

7.2 x 10“ 5 



Elevation 

2.8 x 10~ 4 

0.4 x 10~ 6 

25.2 

1.52 

0.0438 

4.1 x 10“ 4 

2.5 x 10- 5 

220 

32 


to 

to 




to 

to 




28 x 10- 4 

4 x 10~ 6 




11.8 x 10- 4 

7.2 x 10- 5 




5. Hydraulic Line. There are two lines: A and B. A model for line A is developed, and the model 
for line B is similar (index “a” should be replaced with ,s b”). Line A has four inputs, flows q a , q a v,q a pv, 
and q asv , and a single output, pressure p a (refer to Fig. 5). From [7], one obtains the line- A model as an 
integrator, with the negative feedback signs as in Fig. 5: 


Pa — k[ a ( q a -f- (Jav qap i/as) 


(31a) 


Similarly, the line-B model is obtained: 

Pb = kib{ qb qbv — qbp qbs} 


(31b) 


As before, defining p = p a — Pb, one obtains 


p = ki{-q + q v - q p - q s ) 


(32) 


In these equations, the gains are 



(33) 


where /3 is the effective bulk modulus (capacitance of the line), /? = 1.29 x 10 6 kPa (1.87 x 10 5 psi), and 
v Q is the total volume, v Q — 27, 200 cm 3 (1660 in. 3 ), so that ki — 47.3 kPa/cm 3 (113 psi/in. 3 ). The values 
are collected in Table 3. 
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Table 3. Parameters of line A and line B. 


Drive 

/3, kPa 

13, psi 

v 0 , cm 3 

v 0t in. 3 

Line A 

Azimuth 

Elevation 

1.29 x 10 6 
1.29 x 10 6 

1.87 x 10 5 
1.87 x 10 5 

27,200 

24,500 

1653.6 

1498.1 

Line B 

Azimuth 

Elevation 

1.29 x 10 6 
1.29 x 10 6 

1.87 x 10 5 
1.87 x 10 5 

27,700 

24,600 

1690.4 

1499.9 


6. Hydraulic System Model. The model of the hydraulic system is derived by combining its 
elements (active servo valve, shorting valve, passive valve, hydraulic lines, and hydraulic motors). By 
introducing the new differential variables, the block diagram in Fig. 5 is simplified to the one in Fig. 6. 
A detailed block diagram of the hydraulic system is shown in Fig. 7. Combining Eqs. (12), (14), (16), 
(32), and (30) (or, alternatively, using the block diagram in Fig. 7), and defining the new state vector 
x h = [xi,X 2 ,X 3 ] t , with three states, xj = q v ,x 2 = q v , and X3 = p, and defining the input current, i, 
motor rate 6 m , and the single-output motor torque, T 0 , one obtains 


Xh — A/,.7'/, + Bho@m T ^ hi i 
7' 0 = Ch^h Dho@m T Dhi j 


(34a) 


where 



Fig. 6. Simplified block diagram of the hydraulic drive. 
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Fig. 7. Detailed block diagram of the hydraulic system. 



(34b) 


The plots of the magnitudes of the transfer function in azimuth and elevation from i to T a are shown 
in Fig. 8. The plots of the magnitudes of the transfer functions in azimuth and elevation from 9 m to T a 
are shown in Fig. 9. 

C. Gearbox Mode! 

The gearbox model was described in detail in [5], and its block diagram is given in Fig. 10. In this 
diagram, T 0 is the motor torque, 9 P is the antenna angular rate, w m is the motor rate, T is the gearbox 
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Fig. 8. Magnitude of the hydraulic drive transfer function from the solenoid current to 

the output torque. 



Fig. 9. Magnitude of the hydraulic drive transfer function from the motor rate to the 

output torque. 



Fig. 10. Block diagram of the gearbox. 


torque, J m is the motor inertia, k g is the gearbox (output) stiffness, and N is the gearbox ratio. This 
model has two inputs, the motor torque, T c , and the wheel (pinion) angular rate, 8 P , and a single output, 
the gearbox torque, T. 

The equations for this system are as follows: 
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Jm&m — T 0 


N 


(35a) 


Denoting the state variables xi = u> m and x 2 = T, one obtains 

• - ~ X2 i ll. 
Xl " NJ m + J m 


(35b) 


(36a) 


± 2 


kgX 1 

“j\T 


kgOp 


(36b) 


Defining the gearbox state as x g = [xi X 2 ] T , input T e and 0 P , and output T and uj m , 
gearbox state-space representation (A g . B g ,C g ): 

Xg = AgXg + B gl T 0 + B g 2 0 p 

T = CglXg > 

@m = Cg2 x g . 


one obtains the 


(37a) 


where 


Ag 



-1 -| 
NJm 
0 




> 


< 7 * 1=10 1 ] 
< 7*2 =[1 0 ] 


(37b) 


D. Drive Model 

The drive model is obtained by combining the state-space representation of the electronic board, 
Eq. (10); the hydraulic system, Eq. (34); and the gearbox, Eq. (37), according to the block diagram in 
Fig. 1. Defining the drive state vector Xd = [xj ,x^,xJ] T , we obtain the state equations 
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> 


(38a) 


id = A d x,i + Bd r r + BdtOp 

T = CdXd 


where 


A d = 


A b 

BhiCb 


0 

Ah 


BiaCgi 


BhoC'g2 F 


N t 

BhiDbeCgi 


N t 


B d r 


B d t 


B g iD h iCb B g iC h A g + B gl D ho Cg 2 + BalD ^ b2C ^ 


Bbi 

BhiDbi 

L B gl D hi D bl 


0 

0 

L B g2 


C d = [0 0 C gl ] 


(38b) 


The plots of the magnitudes of the transfer function in azimuth and elevation from r to T are shown 
in Fig. 11. The plots of the transfer functions in azimuth and elevation from 9 P to T are shown in Fig. 12. 


III. Structure Model 

The structural model is derived from the finite element model of the antenna structure with free 
rotations with respect to the elevation and azimuth axes. The finite element model consists of the 
diagonal modal mass M m (p x p), diagonal natural frequencies matrix f l(p x p), diagonal modal damping 
matrix Z(p x p), and modal matrix $(m x p),p < m, which consists of p eigenvectors fa (mode shapes), 
i = 1, 




(39) 


Let the finite element model have m degrees of freedom, with s inputs u(t), where uissxl vector, 
and with r outputs y(t), where y is r x 1 vector. If the input matrix is B 0 (m x s), the output matrix 
for displacement is C oq (r x m), and the output matrix for rates is C ov (r x m), then the input-output 
relationship is given by the following second-order differential equation: 


Qm T 22fh) m F Qrn — ltl $ B 0 u 

> 

Vm ~^oq^9m F C ov ^Q m 


(40) 


Define the state variable x as follows: 
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Fig. 1 1 . Magnitude of the drive transfer function from the rate command to 
the output torque. 



Xi 


Qm 

_ x 2 . 


Qm 


(41) 


where q m and q m are modal displacements and rates (such that q = $g m ;g is the actual displacement); 
then Eq. (40) can be presented as a set of first-order equations: 


±i = x 2 


x 2 = - fl 2 xi - 2Z£lx 2 + M^QBoUs > 
Vs i C ov (j)X2 ^ 


or in the following form: 


X s AgXg BgUg 

> 

Vs —C s x 


(42) 


(43a) 
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where 


A s = 


0 

-ft 2 


I 

-2 Zft 


B. 


M-WBo 


C s =[C oq $ C ov $] 


(43b) 


is the sought state-space model in modal coordinates. In our case, u s = [T a T e ], where T a and T e 
are torques at azimuth wheels and elevation pinions, respectively. The structure output consists of the 
elevation and azimuth encoder angles and rates, pinion angles, elevation and cross-elevation pointing 
errors, and other structural variables of interest. Two outputs, 9 pa and 9 pe , the pinion rates in azimuth 
and elevation, are of special interest. Thus, the structural state-space equations are as follows: 


x s — A s x s + B sa T a + B se T e 
9pa = (-'pa X s 

8 p e — ('p a X s 

y = C s x j 


(44) 


The modal data obtained from the finite element model consist of 150 natural frequencies, w*; modes, 
<j>i\ and modal masses, m m j, i = 1, • • • , 150. Additionally, based on the measurements, the modal damping 
is assumed to be 1 percent, i.e., Q = 0.01. Based on this information, the state matrix A s , as in Eq. (43b), 
is determined by introducing the matrix of natural frequencies, ft = diag(uJi), and modal damping, 
Z = diag(£i), i — 1, • • ■ , 150. 

The determination of matrices B s and C s is presented here for the azimuth wheel torque input and 
the azimuth wheel rate output. For the azimuth wheel torque input, consider the azimuth wheel of radius 
r a and the azimuth rail of radius R a - Let nodes n i be located at the contact point of the wheel and the 
rail. The torque applied to the wheel generates the force F a at node n\. The force is tangential to the 
azimuth rail. Assuming a rigid pinion, the force F a applied to the wheel is 



(45) 


This force has x and y components, F ax and Fay [see Fig. 13(a)], such that 


X a 

Fax F a COS Oi a — — COS QIq 

r a 


T 

Fay = F a sin a a = — sin a a 

r a 


> 


(46) 
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Fig. 13. Forces and rates at the azimuth pinion: (a) forces and (b) rates. 

and a a is the angle marked in this figure. Let e x and e y denote the unit vector (all but one component 
are zero, and the nonzero component is equal to one), with the unit component at the location of the x 
and y displacement of node n\ in the finite element model. The input, F, to the finite element model 
is F — F ax e x + F ay e y . Therefore, B 0 follows from the decomposition of F, such that F = B a T a . Prom 
Eq. (46), it follows that 


B 0 = cos a a + — sin a a 

r a r„ 


(47) 


Next, from Eq. (43b), it follows that the nonzero (lower) part of B, after introduction of Eq. (47), is 


M l <b T B a = - M 1 cosa a + M* ' '-' y - sin a„ 


_ 1 $ T e x j $ T e y . 

one I A/f ^ oi' 


M m l — cos a a + — sin a 0 

T a T a 


(48a) 

(48b) 


where <j) x and <j> y are vectors of modal components of x and y displacements at node n i : 


0x — ^ Cx — < t > x2i ' ' ' i 0xl5o] 


T \ 


> 


(49) 


^ e y ~ [tfiyh 4 > y2) * * * > 0yl5o] J 

where <\> xt and 4> yi are x and y displacements of mode i at node ni. Therefore, from Eqs. (43b) and (48b), 
one obtains 


B s 


0 

~M m l — cos a a + M" 1 ^ sin a a 
r a r a 


(50) 


The output matrix derivation is presented here for the wheel rate, 6 pa . The wheel rotation is 
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(51) 



Va 

r a 


where v a is the tangential velocity of the wheel at the contact point [see Fig. 13(b)]. If v x and v y are x 
and y components of v a , and a a is the angle marked in this figure, then 


v a = —v x cos a a + Vy sin a a 


therefore, 



— — cos a a + 
r n . 



Q 


(52) 


(53a) 


and in modal coordinates 


^pa 


$el L . V 

— cos ot a H sin a a q m = 

r a r a / 


4>t <t>y \ 

cosa a + — Sin a a q m 

r a r a / 


(53b) 


Finally, the matrix C s , according to Eqs. (43b) and (53b), is 


C s 


4>l 

0 — — cos a a 

TV 


sma„ 


(54) 


The structural model consists of m = 150 modes or 300 states. Modes not participating in system 
dynamics are eliminated. Observability and controllability properties in the balanced representation are 
used to determine insi g nifi cant, modes. The balanced representation [9] is a state-space representation 
with equally controllable and observable states. The Hankel singular value is a measure of the joint 
controllability and observability of each balanced state variable. The stated with small Hankel singular 
values are deleted as weakly excited and weakly observed, causing minimal modeling error. 

For flexible structures with small damping and distinct poles, the modal representation is almost 
balanced, c.f. [10-12], and each mode is considered for the reduction separately. For a structure with 
m modes, matrix B s has 2 m rows, and C s has 2m columns. Denote b s as the last m rows of B s , c q as 
the first m columns of C s , and cv as the last m columns of C s . Then b si is the ith row of b s ,c q i is the 
zth column of c q , and tv* is the ith column of cv. Denote /3^ = & si 6T,a: gi = <? qi c q *, and a ri — c^c ri . The 
Hankel singular value for the ith mode is given in [11] and [12]: 


it 


WbiPsi y/rfiP%i + w li^ 2 i a li 


40? 


(55) 


where the weighting factors wu > 0, w qi > 0, w r i > 0, and i = 1, ■ • • , m. 

Care should be taken when determining Hankel singular values. Units should be consistent; other- 
wise, some inputs or outputs receive more weight in Hankel singular-value determination than necessary. 
Consider, for example, the azimuth encoder reading in arcseconds and the elevation encoder reading in 
degrees. For the same angle, the numerical reading of the azimuth encoder is 3600 larger than the eleva- 
tion encoder reading; hence, the elements for the azimuth output are much larger than those for elevation. 
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On the other hand, some variables need more attention than others: Pointing error and encoder readings 
are the most important factors in the antenna performance; hence, their importance has to be emphasized 
in mode evaluation. For consistency of units and importance of variables, the weighting factors Wbi,w qi , 
and w r i are introduced. Typically, weights are set to 1. 

For each mode, the Hankel singular value is determined and used to decide on the number of modes in 
the reduced structural model. For the rigid body modes, Hankel singular values tend to infini ty; hence, 
rigid body modes are always included in the reduced model. Hankel singular values of the 150 modes of 
the antenna model are plotted in Fig. 14. The reduced order model consists of 24 modes: 2 rigid-body 
modes and 22 flexible modes. 

The plots of the transfer function in azimuth and elevation (magnitude and phase) from the wheel 
(pinion) torque T to the axis rate 6 are shown in Fig. 15. They show that the azimuth transfer function 
has low frequency resonances (about 1.2 and 2.2 Hz), which are absent in the elevation transfer function. 




Fig. 15. Magnitude of the transfer function of the antenna structure: (a) direct coupling 

and (b) cross-coupling. 
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SV, Rate Loop Model 

A rate loop block diagram is presented in Fig. 16, where T a and T e denote the drive torques, 0 pa and 
6 pe denote pinion rates, and r a and r e are rate commands in azimuth and elevation, respectively. The 
state-space equations are combined from the state equations of the azimuth and elevation drives [see 
Eq. (38) and add subscript “a” for the azimuth drive and subscript “e” for the elevation drive] and the 
structure [see Eq. (44)]. Combining them, and defining the rate loop state vector x r as x r = [xd a , x de ,x s }, 
where Xda and Xd e are azimuth and elevation drive states, one obtains the rate-loop state-space equations: 


Xr = A r X r + B ra'T a + B 

re T e 


y — C r x r 

► 

0 a = CdXf 

Q e = C e x r j 


(56a) 


where 


- -uu 


A r — | 0 Ad e Bdte^pe 

_ B sa Cda B se C de A s 


B r a — 


Bdra 

0 

0 


B r 


0 

B d re 

0 


Cr = [0 0 C s 


(56b) 


where 6 a and 9 e are azimuth and elevation encoder readings, C pa and C pe are the output matrices for the 
azimuth and elevation pinion rates, and C a and C e are the output matrices for the azimuth and elevation 
encoders, respectively. 

Figure 17 shows the magnitude of the transfer function from the azimuth rate input r a to the azimuth 
encoder rate 9 a (solid line) and the magnitude of the transfer function from the elevation rate input r e to 
the elevation encoder rate 9 e (dashed line). The figure shows that the required identity relationship for 
low frequencies is not acquired. The magnitude of the transfer functions for frequencies less than 0.3 Hz 
is 0.74, below the required 1, due to inaccuracy in the model parameters (mainly in the hydraulic part). 
T his drawback can be removed by the experimental investigation of the parameters of the hydraulic 
drives, such as motors, valves, and lines. However, this inaccuracy is corrected by the position feedback 
loop, as will be shown later. The high-frequency peaks in azimuth and elevation (8 Hz in azimuth and 
20 Hz in elevation) are the gearbox resonances. 
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Fig. 16. Block diagram of the rate loop model. 



Fig. 17. Magnitude of the transfer function of the rate loop model. 


V. Position Loop Model 

The rate loop system with the proportional and integral (PI) controller is shown in Fig. 18, where e a 
and e e are the azimuth and elevation servo errors. For the series connection of the rate loop system and 
the controller, as in Fig. 18(a), define the state vector = [x ai x ei xj t ] with the new state variables 
x ei and x a i (integrals of the errors) such that 


ai — 1 



(57) 


The system output y is defined in Eq. (56a), the encoder output is 0 T = [6 a 0 e \, and the input is 

e T = [e a e e ]. The inputs to the rate loop systems are obtained from Fig. 18(a): 


X a — k pa &a T ^ 2 a X a i ) 


V e — k pe e e “I - ki e X e i 




(58) 


where k pe ,k ie ,k pa , and ki a are proportional and integral parameters of the controllers. Combining the 
equations for the rate loop system with Eqs. (57) and (58), one obtains 
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(a) 




Fig. 18. Position loop: (a) open and (b) closed. 


where 


Xq — A. 0 Xq “I - B 0 c 

9 = C 0 x 0 > 

y = Cx 0 


0 

0 

kia Bra 


0 

0 

kieB r e 


0 

0 

At 


Bo = 


I 

0 


0 

7 


kpaBra kp e I3re _ 


Co 


0 0 c a 
0 0 Ce 


C=[ 0 0 C r ] 


For the closed-loop system [see Fig. 18(b)], 


(59a) 


(59b) 


e = c — 6 (60) 

where c T = [c a c e ] is a command signal in azimuth, c a , and in elevation, c e . Introducing Eq. (60) to 
Eq. (59), one obtains 


x c i — x qi -f- B r c 


y = Cx c i 


(61a) 


where 
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A-ci — A 0 B 0 C c 


(61b) 


The simulations shown in Figs. 19 through 22 have two sets of assumptions: (1) a proportional gain 
of 1 in azimuth and elevation and an integral gain of 0.3 and (2) a proportional gain of 0.7 in azimuth 
and elevation and an integral gain of 0.2. The closed-loop transfer functions from azimuth command 
to azimuth encoder are shown in Fig. 19(a), and those from elevation command to elevation encoder 
are shown in Fig. 19(b). They show a bandwidth of 0.1 Hz. The cross-coupling transfer functions from 
azimuth command to elevation encoder and from elevation command to azimuth encoder are shown in 
Fig. 20. They show low-level cross-coupling. The closed-loop step responses from azimuth command 
to azimuth encoder are shown in Fig. 21(a), and those from elevation command to elevation encoder 
are shown in Fig. 21(b). They show a 20- to 30-percent overshoot and a 7- to 9-s settling time. The 
cross-coupling from azimuth step command to elevation encoder and from elevation step command to 
azimuth encoder is shown in Fig. 22. The cross-coupling is of the order 10 -3 . 
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Fig. 1 9. Magnitude of the transfer function of the 
closed-loop system: (a) azimuth command to 

azimuth encoder and (b) elevation command to 
elevation encoder. 
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Fig. 20. Magnitude of the transfer function of the 
closed-loop system: azimuth command to 

elevation encoder and elevation command to 
azimuth encoder. 



01 23456789 10 

TIME, s 


Fig. 21. Step responses of the closed-loop 
system: (a) azimuth command to azimuth 

encoder and (b) elevation command to elevation 
encoder. 





Fig. 22. Step responses of the closed-loop 
system: azimuth command to elevation encoder 
and elevation command to azimuth encoder. 


VI. Wind Disturbance Simulations 

Wind gust disturbances were modeled similarly to the DSS-13 antenna (see [13]) using the wind tunnel 
pressure distribution on the dish taken from Blaylock. 1 Their time history is generated using the wind 
Davenport spectrum (see [14] and [15]), determined for the Goldstone site. The simulations for the 
50 km/h wind gave the results listed in Table 4 and compared with the simulation results of the DSS-13 
antenna. The table shows that DSS 14 has better disturbance rejection properties (at the encoders) than 
has the DSS-13 antenna. 


Table 4. Servo errors in mdeg (3cr rms) 
for 50 km/h wind gusts. 


Drive 

Front wind 

Side wind 

Elevation, DSS 14 

2.6 

0.7 

Elevation, DSS 13 

14.6 

1.9 

Azimuth, DSS 14 

0.1 

2.1 

Azimuth, DSS 13 

0.5 

2.3 


VII. Conclusions 

An analytical model of the DSS-14 antenna has been developed. The rate loop model consists of the 
structural model (derived from the finite element model), gearbox model, hydraulic servo, and electronic 
boxes. The position loop was closed, and the time and frequency responses were simulated. The wind 
pointing errors of the DSS-14 antenna have been simulated. The model allows for detailed simulation of 
antenna dynamics and for modifications and improvements to the antenna control system. 

The simulations confirmed that the use of encoders located at drives limits the performance of the 
antenna (mainly by reducing its bandwidth to 0.1 Hz). The use of the master equatorial or new encoders 


1 R. B. Blaylock, “Aerodynamic Coefficients for Model of a Paraboloidal Reflector Directional Antenna Proposed for a 
JPL Advanced Antenna System,” JPL Interoffice Memorandum CP-6 (internal document), Jet Propulsion Laboratory, 
Pasadena, California, 1964. 
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located close to the axes of rotation of the antenna (similarly to the 34-m antennas) would allow expansion 
of the bandwidth to 0. 7-1.0 Hz. 


The antenna model needs further improvement. First, in this model, certain parameters of the hy- 
draulic drive are known with rather poor accuracy, and it influences the accuracy of the antenna model. 
It is essential to use experimental techniques to get more precise values of the parameters. Secondly, the 
RF pointing errors (in elevation and cross-elevation) of the antenna should be determined in order to 
evaluate the precision of the antenna pointing. 
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Appendix 

Transfer Function Derivation 


Each component of the electronics board is composed of operational amplifiers (opamps), resistors, 
and capacitors. The basic configuration of an inverting opamp circuit is shown in Fig. A-l. The “+” 
terminal of the opamp is grounded; thus, the “— ” terminal voltage is zero, called a virtual ground. In 
this situation, the currents ii and i 2 flowing through impedances Z\ and Z 2 are equal to 


} (A-l) 


and their sum is zero; that is, ii = — i 2 . Introducing them to Eq. (A-l) gives 






*2 


Vout 

z 2 


Vout — hv^n 


(A-2a) 


where 



(A-2b) 
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Fig. A-1 . Opamp circuit. 


I. Transfer Function Gf 0 

A schematic for the transfer function G to is shown in Fig. A-2(a), where the notation and the value 
of each element were taken from JPL Drawing 9479871D. 2 This schematic can be simplified to the one 
shown in Fig. A-2(b). In this figure, 

R s = {^R^q T T + R 59 ^) + Fi ’)2 — 49.5 kQ (A-3a) 

where R 5e = i ? 57 = R58 = R59 = 100 kSl, and Rq2 = 24.5 kfi; thus, R s = 49.5 kfi. The component Z\ is 

^1 = i?63 + = ^3 + i?64 - 91.1 m (A-3b) 

i + rt64^40 s 

where R$s = 40 kO, R $ 4 = 51.1 kf2, and C 40 = 0.15 fiF. The time constant R 64 C 40 — 0.0077 s is small, 
thus neglected. Denote Rmta — 9.7 kfl and R rn t e = 7.8 kfi the motor resistances in azimuth and elevation, 
respectively, and C 41 = 0.15 /jF. Then, the component Z 2 for the azimuth drive is as follows: 

Z 2 = 7 7 p rr,t V - - R mta = 9.7 kfi (A-4a) 


and for the elevation drive, 


Z 2 = i T ~ ^ = 7 ‘ 8 (A-4b) 

1 -r -fL 7 n f :e 04 xS 

The time constants R m taC 4 \ = 0.0015 s and R m teC 41 = 0.0012 s are of the order 10~ 3 s, thus considered 
small, and neglected. 

The transfer function G t0 for azimuth is 


Gto — 


R P 

R s + Rp 


8800 

49, 500 + 8800 


0.151 


(A-5a) 


and for elevation, it is 


2 JPL Drawing 9479871D (internal document), Jet Propulsion Laboratory, Pasadena, California. 
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(A-5b) 



Fig. A-2. Schematic for the transfer function G/q: (a) full and (b) simplified. 


II. Transfer Functions G r i and G r 2 

The transfer functions G r \ and G r 2 are determined simultaneously. Their schematic is given in 
Fig. A-3(a), the parameters parameters of which are 

R i5 = 750 kft 

R 50 — 100 kfi 

i?5i = 12.1 kfl 

R 52 = 442 kfl 

R 53 = 442 kO 

Rq5 — 909 kO 

i?66 = 90.9 kfl 

C 31 — 1 fJ-F 

C 42 = 0.1 fiF 

The schematic from Fig. A-3(a) can be transformed to the form shown in Fig. A-3(b). The value of Z 3 
is as follows: 

Zz = i?66 + 1 p 65 /^ — -^66 + #65 = 10 6 kfl (A-7) 

1 + i<65W2S V 
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In this variable, the small time constant ResC^ = 0.0909 s was ignored. 


The value Z 4 is obtained as 


Z 4 = R 53 4- 


Rp 

1 + R 0 C 31 s 


- 4.65 x 10 6 


1 + 0.400s 
1 + 4.205s 


(A-8) 


where R 0 = (R50 + R52 + R50R52)/ R51 = 4205 kO. 

Having determined Z 3 and Z 4 , the transfer functions G r 1 (from v r to v s ) and G r 2 (from v to to v s ) are 
obtained: 


G r i=~ = 6.20 G a 

“15 

G T 2 = — ~ = —4.65 G 0 
63 


> 


(A- 9 ) 


where 


1 + 0.400s 
0 ~ 1 + 4.205s 


(A-10) 


is the transfer function of a lag compensator. 




Fig. A-3. Schematic for the transfer functions and G^: (a) full and (b) simplified. 


III. Transfer Function G s 

The transfer function G a is determined from the schematic in Fig. A-4(a), and is shown in compact 
form in Fig. A-4(b). For this schematic, Ri 3 = 100 kO, R 36 = 10 kQ, R43 = 24.9 kfl, and Gis = 0.1 /iF; 
therefore, one obtains 


^5 


R36 

1 + R36G18S 


“ R 36 = 10 kO 


(A-ll) 


where i2 36 Ci8 = 0.0015 s = 0. Since Vi — v s Z 5 /R 43 , and i s = i>i (Z 5 1 + R 1 ), thus, 
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A Prototype Ka-/Ka-Band Dichroic Plate With 
Stepped Rectangular Apertures 

J. C. Chen, P. H. Stanton, and H. F. Reilly, Jr. 

Communications Ground Systems Section 


A prototype five-layer Ka- /Ka-band dichroic plate was fabricated and measured. 
This dichroic plate was designed to pass Ka-band uplink (34.2-34.7 GHz) and to 
reflect Ka-band downlink (31.8-32.3 GHz) for dual-frequency operation in the Deep 
Space Network to support the future Cassini mission. The theoretical calculation 
and the experimental measurement of the reflected resonant frequencies were within 
0.24 percent for circular polarization. The computer program, which was used to 
design the dichroic plate with stepped apertures, was then verified. 


I. Introduction 

A dichroic plate was needed for simultaneous Ka-band uplink and downlink operation. In order to 
diplex two frequency bands with only a 1:1.07 ratio, stepped apertures were chosen for the Ka-/Ka-band 
dichroic plate design [1], The stepped apertures, acting as resonator filters, pass the high-frequency band 
(Ka-band uplink, 34.2-34.7 GHz) and reflect the low-frequency band (Ka-band downlink, 31.8-32.3 GHz). 
The five-step aperture is a rectangular waveguide with two thin irises that divide the waveguide into three 
sections (Fig. 1). 


II. Calculated Reflected and Transmitted Radiation Patterns and Power 

The reflected and transmitted radiation patterns and power were calculated based on an incident horn 
pattern model at 32.0 and 34.5 GHz [2]. A 26-dB horn radiation pattern was used as the input pattern 
to the Ka-/Ka-band dichroic plate. The reflected and transmitted patterns were calculated at every 
4> = 12.25-deg interval and 9 = 1-deg increment. The 9 is the angle measured from the axis of rotation of 
the feedhorn (or its image in the case of a reflected pattern). The <p is the rotational angle around the above 
axis. The Cray supercomputer was used to perform this intensive calculation since 80 waveguide modes 
were required to ensure accuracy. The reflected and transmitted patterns are shown in Figs. 2 through 5 
for the two linear polarizations that are orthogonal to each other at 32.0 and 34.5 GHz, respectively. The 
transmitted and reflected patterns of 32.0 GHz are similar to the incident horn pattern except that the 
peaks of the transmitted patterns are about —45 dB. The transmitted patterns of 34.5 GHz are similar 
to the incident horn pattern, while the peaks of the reflected patterns are about —20 dB. 

The transmitted and reflected power can be calculated by integrating the transmitted and reflected 
patterns. It was found that 99.99 percent of the power is reflected and 0.01 percent leaks through the 


143 



144 


y 



Fig. 1. Geometry c 
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plate at 32 GHz. The power leakage may contribute 0.03 K of the noise temperature to the antenna 
system at DSS 13 at 34.5 GHz; 2.0 percent of the power is reflected and 98 percent is transmitted. The 
reflection and transmission coefficients versus frequencies were also calculated. The detailed information 
is in [1]. 
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Fig. 2. Calculated reflected and transmitted radiation pattern of a 
Ka-/Ka-band dichroic plate at 32 GHz for linear polarization. 
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Fig. 3. Calculated reflected and transmitted radiation pattern of a 
Ka-/Ka-band dichroic plate at 32 GHz for the orthogonal linear 
polarization. 
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Fig. 4. Calculated reflected and transmitted radiation pattern of a 
Ka-/Ka-band dichroic plate at 34.5 GHz for linear polarization. 



0, deg 


Fig. 5. Calculated reflected and transmitted radiation pattern of a 
Ka-/Ka-band dichroic plate at 34.5 GHz for the orthogonal linear 
polarization. 


III. Fabrication 

The prototype Ka-/Ka-band dichroic plate was a 177.8-mm-by-177.8-mm five-layer copper plate with 
a 152.4-mm-diameter perforated area (Fig. 6). The advantage of having only two different aperture sizes 
is a reduction of the fabrication cost. Since multiple metal sheets can be stacked together to be wire 
electrical-discharge machined when the sheets have an identical pattern, only two sets of sheets need to 
be run though the machine. 


Previous dichroic plate fabrication dictated that, due to the tight tolerance and the cost, the wire 
electrical-discharge machine was the best way to fabricate rectangular apertures. The five-layer plate 
required a way to combine the layers into a solid structure. Holes in all four corners of the five plates 
were used to align the rectangular apertures and the five plates. 
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Fig. 6. Prototype Ka-/Ka-band dichroic plate. 


Pure silver was deposited on both sides of the two thin plates in an electroless solution for a thickness 
of from 0.00508 to 0.00762 mm. The plates were then aligned with the four aligning pins originally used 
to locate the rectangular apertures. The fusing was accomplished in a reduced atmosphere of a hydrogen 
oven. A test sample made from regular copper was not adequate for this technique since the oxygen in 
the copper went to the surface during the brazing process. This test result dictated that only oxygen-free 
copper (less than 0.001 percent) was acceptable. 

After visual inspection, the five-layer plate was mechanically inspected using the Monarch, Cortland, 
and Bafire Probe technique. Data points of the location and size of the rectangular apertures were 
recorded and averaged. The inspected dimensions of the apertures turned out to be 3.8331 mm by 
3.9406 mm (design dimensions of 3.8100 mm by 3.9268 mm) for small apertures and 4.6200 mm by 
4.6116 mm (design dimensions of 4.6126 mm by 4.6228 mm) for the large apertures. All dimensions were 
within a 0.0254-mm tolerance. 


IV. Measurement 

The reflection coefficients of the dichroic plate were measured on a Hewlett Packard 8510 network 
analyzer (Fig. 7). Two orthogonal linear polarizations, TE and TM, were measured. (When the electrical 
field of the wave is perpendicular to the plane of incidence, it is called TE polarized; if the magnetic 
field of the wave is perpendicular to the plane of incidence, it is called TM polarized.) A 22-dB cor- 
rugated horn was located 104-mm away from the dichroic plate along a 30-deg angle of incidence. A 
second horn was placed 1073 mm from the dichroic plate along the direction of the reflected beam. The 
reflection response curve shows two resonant frequencies for TE and TM linear polarizations, respec- 
tively (Figs. 8 and 9). The theoretical values were recalculated for the actual dimensions of the proto- 
type dichroic plate. Results for circular polarization were computed by using the two orthogonal linear 
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Fig. 7. Experimental setup for the reflection measurement. 
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Fig. 8. Measured and calculated reflection coefficient versus 
frequency for TE linear polarization. 


polarizations (TE and TM) and the phase difference between them (Fig. 10). The calculated and measured 
resonant frequencies are less than 0.25 percent in error for circular polarization (Table 1). 


A transmission measurement for TE and TM polarizations was also made in order to check the perfor- 
mance of the reflecting band (31.8-32.3 GHz) (Fig. 11). The 22-dB transmit horn was located 127.0 mm 
away from the dichroic plate, which was tilted 30 deg from normal. The five-layer dichroic plate is ap- 
proximately three times thicker than a one- layer dichroic plate such as the X-/Ka-band dichroic plate 
[3]. Therefore, the transmitted beam path was not the same as the incident beam path. Consequently, 
the receiving horn was offset 13.97 mm from the incident beam path. For reference, the transmission 
coefficients were first measured without the dichroic plate. Then the transmission coefficients were mea- 
sured again with the dichroic plate between the two horns. The measured and calculated transmission 
coefficients showed good agreement for TE and TM polarizations, respectively (Figs. 12 and 13). The 
transmitted power was below —35 dB at the Ka-band downlink for circular polarization (Fig. 14). 




V. Conclusion 

The prototype five-layer Ka-/Ka-band dichroic plate was successfully fabricated. The fabrication 
technique was proven adequate for the multilayer dichroic plate. The experimental results and theoret- 
ical predictions showed good agreement for the resonant frequency, within 0.25 percent. Furthermore, 
the computer code for analyzing the dichroic plate with stepped rectangular apertures was verified. A 
Ka- /Ka-band dichroic plate will be fabricated and implemented for an X-/X- to Ka-/Ka-band demon- 
stration at DSS 13. 
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Table 1. The measured and calculated reflection resonant frequencies 
of the prototype Ka-/Ka-band dichroic plate. 


Polarization 

Resonant 

frequency 

Measurement, 

GHz 

Calculation, 

GHz 

Error 

percentage 

TE linear polarization 

1st 

33.56 

33.625 

0.19 


2nd 

34.37 

34.5 

0.38 

TM linear polarization 

1st 

33.41 

33.485 

0.22 


2nd 

34.32 

34.40 

0.23 

Circular polarization 

1st 

33.47 

33.55 

0.24 


2nd 

34.36 

34.44 

0.23 


127.00 mm 



EDGE DIFFRACTION 


Fig. 11. Experimental setup for the transmission measurement. 



Fig. 12. Measured and calculated transmission coefficient 
versus frequency for TE linear polarization. 
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In February 1995, seq uence-of- even ts (SOE)-driven automation technology was 
demonstrated for a Voyager telemetry downlink track at DSS 13. This demonstra- 
tion entailed automated generation of an operations procedure (in the form of a 
temporal dependency network ) tom project SOE information using artificial intel- 
ligence planning technology and automated execution of the temporal dependency 
network using the link monitor and control operator assistant system. This article 
describes the overall approach to SOE-driven automation that was demonstrated, 
identifies gaps in SOE definitions and project profiles that hamper automation, and 
provides detailed measurements of the knowledge engineering effort required for 
automation. 


I. Introduction 

The Voyager 1 spacecraft is cruising at 17.5 km/s toward the outer edge of the solar system. Though 
its onboard systems are mostly asleep during this phase of its mission, Voyager’s health metrics are 
continually sent to Earth via a telemetry signal radiated by its 40-W transmitter. It will take 8 hours 
at the speed of light for the signal to reach its destination, Earth, a billion miles away. Upon arrival, 
the telemetry signal is received by an extremely sensitive ground communications system, NASA’s Deep 
Space Network (DSN), where it is recorded, processed, and sent to the mission operations and Voyager 
project engineers, who assess the health of the spacecraft based on the contents of the signal. 

The type of activity just described occurs daily for dozens of different spacecraft and projects. Though 
the process of sending signals from a spacecraft to Earth is conceptually simple, in reality there are many 
Earth-side challenges, both technical and institutional, that must be addressed before a spacecraft’s signal 
is acquired and transformed into useful information. The purpose of this article is to identify some of these 
challenges and propose ways of automating the process of capturing data from a spacecraft. In particular, 
this article focuses on how to automatically transform a flight project service request into an executable 
set of DSN operations that fulfill the request. To understand in more concrete terms the processes 
we propose to improve, consider first the following example of the typical activities that currently occur 
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Fig. 1. Partial process map for network preparation and data capture processes. 

when the DSN provides a service in response to a project request. Figure 1 shows a partial process map 
of the steps taken to fulfill the service request. 1 

A team within the Voyager flight project decides to perform a periodic check of the health status of 
the Voyager 1 spacecraft. Since the spacecraft’s health metrics are sent via a telemetry signal, a project 
representative submits a service request to the Telecommunications and Mission Operation Directorate 
(TMOD) (see Fig. 1) for a telemetry downlink activity to be performed, and the time, day, and equipment 
are negotiated with a representative from the DSN. Once the activity has been scheduled (Fig. 1, “schedule 
resources”), the Voyager team develops a project sequence of events (PSOE), 2 which is a time-ordered 
description of the events that should take place during the activity (Fig. 1, “generate PSOE”). The 
PSOE is written in a high-level language 3 that provides a set of keywords and parameters to specify 
DSN functions and spacecraft events. Some of the keywords used in the SOE are actions that must 
be performed by the DSN, e.g., configure the receiver or acquire a downlink. Other keywords describe 
actions that will be taken by the spacecraft, e.g., a change in bit rate or the modulation index during the 
track activity. Once completed, the Voyager team sends the SOE to DSN operations, whose job it is to 
implement the requested activities. 

After DSN operations receives the SOE keyword file, which is a high-level description of the events and 
activities that must be performed, operations personnel supplement the PSOE with additional information 


1 This is a modified version of the full process maps found in Final Report of the Reengineering Team (RET) for Services 
Fulfillment (Operations) (internal document), Telecommunications and Missions Operations Directorate (TMOD), Jet 
Propulsion Laboratory, Pasadena, California, March 14, 1995. 

2 This is a simplified description of the actual process. In reality, a project develops an integrated PSOE describing all of 
the activities planned for its spacecraft for a 2-week period. The integrated SOE is usually developed using the project’s 
own language and later converted into the DSN keyword format described above. There are currently some efforts within 
the DSN to standardize the set of services it offers to its customers. With all of the translation required to go from the 
project language to the DSN keyword format, it seems desirable to develop an SOE language or tool that can be reused 
by multiple flight projects and that also would eliminate the need to translate into the DSN format. This is currently 
being addressed by the SOE reengineering team. 

3 Flight Project Interface to the DSN for Sequence of Events Generation, DSN Document 820-13, OPS-6-13 (internal 
document), Jet Propulsion Laboratory, Pasadena, California, August 31, 1995. 
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concerning resource allocation and activity-specific configuration (Fig. 1, “generate TSOE” and “generate 
data types and formats”). 4 The TSOE (which stands for Telecommunications and Mission Operations 
Directorate Ground Network sequence of events) is sent to the deep-space communications complex, 
where operators read and transform its contents into a set of operations (ops) procedures, which, when 
executed, will satisfy the goals stated in the TSOE (Fig. 1, “connection ops”). Based on the keywords 
in the TSOE, the operator chooses procedures to accomplish the goals of the activity and executes them 
by manually typing and sending command directives to the subsystems in the assigned communications 
link. The execution of the issued directives is monitored by the operator, who closes the control loop by 
determining whether the directive was successful or not. By “successful” execution, we mean that the 
directive had its intended effect on the subsystem. Thus, to successfully transform the TSOE from a list 
of time-tagged keywords into an operational procedure, the following steps are performed: 


(1) Determine the appropriate directive(s) and parameters to achieve the current keyword 
goal 

(2) Determine that the necessary conditions required for the execution of the directive are 
satisfied (i.e., verify preconditions of the directive) 

(3) Execute the directive (manual and computer-initiated activities) 

(4) Verify that the directive was correctly received by the subsystem 

(5) Verify that the desired subsystem state has changed (verify successful accomplishment 
of the postconditions of the directive) 

These are the steps taken under normal cheumstances. In addition, the operator monitors the overall 
progress of the activity and the health status of the subsystems by viewing a host of different computer 
displays. Exceptions must be handled in real time by the operator, who invokes recovery procedures when 
there are equipment anomalies or failures. Thus, from the start to the finish, the process of capturing 
data from a spacecraft currently involves many manual steps. 5 


II. Automated DSN Operations — The Concept 

While current DSN operations are extremely labor- and knowledge-intensive, the DSN Technology 
Program has been investing in the groundwork for a new model of operations that automates major 
portions of the DSN data capture and network preparation processes. In this new model, generic services 
are electronically requested of the DSN by a project. These services are then provided to the project in 
a manner completely transparent to the project. Within the DSN services fulfillment process, routine 
operations will be completely automated, greatly enhancing the flexibility and responsiveness of the DSN. 
Furthermore, many of the costly, knowledge- and labor-intensive processes will be automated, resulting 
in a more efficient, streamlined DSN operations structure. 

The Network Automation Work Area, funded by the DSN Technology Program, has been develop- 
ing and demonstrating technology to fulfill this vision of automated data capture and network prepa- 
ration. This work and these demonstrations can be described in terms of three areas (see Fig. 2): 
resource allocation — assignment of global DSN resources necessary to fulfill the requests; procedure 


4 Again, this is a simplified account of the actual process; the project SOE keyword file requires additional information 
before it is sent to the deep-space tracking station for implementation. The SOE also triggers the generation of predict 
and support data needed by various subsystems during the activity, e.g., the antenna needs pointing predicts so that it 
will be oriented to receive the spacecraft signal. 

5 The complete process map is in DSN Document 820-13, OPS-6-13, op cit. 
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Fig. 2. Automating the DSN's data capture and network preparation processes. 


generation — generation of appropriate executable procedures to fulfill specific service requests in light of 
the allocated equipment; and procedure execution — executing the appropriate procedure and ensuring 
that the appropriate services are correctly fulfilled. Below we outline in further detail ongoing efforts in 
each of these areas. 

In the resource allocation area, the OMP-26 scheduling system has been fielded to support scheduling 
of the 26-m DSN antenna subnetwork. Current efforts involve extending the OMP-26 system to schedule 
further classes of DSN assets (extension beyond the 26-m subnet); expansion of OMP-26 coverage to the 
real-time area, where changing requests and equipment availability force modifications of the operational 
schedule; and further extension and maturation of the OMP-26 automated scheduling capability. 

In the procedure generation area, we have been adapting planning technology that has been proven 
and fielded in science data processing applications to the problem of automated generation of antenna 
operations procedures. This work involves development of the DSN planner (DPLAN) system, which 
generates an executable antenna operations procedure, called a temporal dependency network (TDN) [3], 
from a high-level track specification and equipment assignment. 

In the procedure execution area, the link monitor control operator assistant (LMCOA) system is being 
further refined and generalized to support automated tracking using DSN antennas and subsystems. The 
LMCOA system provides closed-loop control, execution monitoring, and error recovery to provide robust 
procedure execution for data capture. LMCOA has been used for the Ka-Band Antenna Performance 
(KaAP) experiments, the Jupiter Patrol, and the phased-array experiment at DSS 13. 
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III. SOE-Driven Automation Demonstration 

A. Overview and Objectives 

In February 1995, a comprehensive demonstration was conducted to validate the concept of end-to- 
end SOE-driven automation. In this demonstration, the procedure generation and procedure execution 
elements of automation were demonstrated by performing a Voyager telemetry downlink track using 
the 34-m beam-waveguide (BWG) antenna at DSS 13 in Goldstone, California. There were four major 
objectives to this demonstration: 


(1) To demonstrate that an executable operations procedure could be generated using the 
DPLAN automated planning system [7]. 

(2) To demonstrate that the TDN generated by DPLAN could be correctly executed by the 
LMCOA procedure execution system to successfully perform the requested data capture 
process. 

(3) To further understand the institutional and infrastructure difficulties in fielding tech- 
nologies for (1) and (2). Specifically, what types of information are currently missing 
from SOEs and affiliated request information that would be needed to automatically 
generate, parameterize, and execute TDNs? Additionally, how much effort in knowledge 
engineering was required to endow DPLAN and LMCOA with the capability for this 
track? How much effort would it take to extend to a related track for a different project? 
To a different type of track? 

(4) To acquire both qualitative and quantitative data on the amount of effort required to 
perform knowledge engineering to field systems for tasks (1) and (2) above [4,6,8,10]. 

The remainder of this section outlines at a high level how each of these capabilities was demonstrated 
and describes the interrelations between different elements of the demonstration. The remaining sections 
in this article describe the underlying technologies and representations. 

The first goal of the end-to-end demonstration was to show automated generation of DSN operations 
procedures. In this element of the demonstration, DPLAN used the PSOE, equipment allocation, and a 
project profile file to construct a TDN that could be executed to perform the track. Generating the TDN 
involves 


(1) Determination of the correct TDN blocks (to achieve the tracking goals in the context 
of the particular equipment assignment and project profile) 

(2) Determination of the correct dependencies among blocks (e.g., block A must be completed 
before block B) 

(3) Determination of the correct parameters for each block (in the context of the tracking 
goals, equipment assignment, and project profile) 

Some elements of TDN parameterization are deferred until plan execution time (to be performed by 
LMCOA). DPLAN uses information about how to achieve track goals using generic TDN blocks in the 
context of various equipment configurations and project parameters. This modular representation of op- 
erations procedures allows for a more compact knowledge base, removing the need to maintain a complete 
set of complex end-to-end TDNs to cover all combinations of track types, equipment configurations, and 
project classes — a set that may number in the hundreds. Each TDN block is a procedure executable by 
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the procedure execution engine, LMCOA, and consists of a sequence of subsystem directives tailored to 
achieve a higher-level goal in context of a specific type of track. 

The second goal of the end-to-end demonstration was to exercise the automated procedure execution 
capability of the LMCOA system on the TDN produced by the DPLAN. This involved demonstrating 
the data capture capability of the LMCOA for a series of Voyager 2 downlink telemetry passes at DSS 13. 
While LMCOA has been in use since the fall of 1994 for KaAP and Jupiter Patrol tracking activities, 
tracking Voyager 2 represented the first time LMCOA had been used to track a spacecraft. Thus, one 
point of the demonstration was to show the generality of the LMCOA approach to both spacecraft and 
nonspacecraft tracking activities. 

The third goal of the end-to-end demonstration was to characterize and quantify the effort required 
to implement the automation architecture. This effort can be attributed to two factors. The first part 
of the implementation effort involved determining elements of information and infrastructure that would 
impede implementation of an end-to-end automation system as demonstrated. This goal was furthered 
by requiring all types of project, spacecraft state, and sequence-related information to be contained in 
the project SOE and project profile. Because the project SOE is fixed, all mission information would 
need to be contained in the project profile. Thus, examining the types of information in the project 
profile provided valuable insight into information missing from the project SOEs. The second part of 
the effort was to determine how much knowledge engineering was required to implement the automation 
architecture. A careful accounting was made of the knowledge engineering effort required to construct 
and validate the DPLAN knowledge base (used in generating the TDN), the TDN blocks themselves, the 
parameterization of the TDN blocks by both DPLAN and LMCOA, and the directives that make up the 
TDN blocks used in the demonstration. 

We think that the labor estimates we derived in constructing the TDNs represent a realistic estimate 
of the true effort required to construct operation TDNs. The labor estimates reflect generation, manual 
evaluation, and operation validation of the TDNs. All three of these steps are required to produce a 
TDN capable of providing reliable operations support. Other TDNs that have been constructed by other 
groups, while carefully constructed and manually verified, have not yet been validated by usage in an 
automation system. 

B. Demonstration Context: Deep Space Station 13 

Located at Goldstone, California, Deep Space Station 13 (DSS 13) is the deep-space tracking station 
dedicated to research and development activities of the Deep Space Network. DSS 13 has functional 
elements (see Fig. 3) comparable to DSN tracking stations: a 34-m BWG antenna, which is controlled 
by the antenna pointing subsystem; feed horns for S-band (2.3-GHz), X-band (8.45-GHz), and Ka-band 
(32-GHz) signals; a total power radiometer (TPR) for measuring broadband signal strength; a weather 
subsystem for tracking wind, humidity, and other atmospheric conditions affecting signal acquisition; a 
microwave switch controller subsystem; and a receiver, the TP-13 telemetry processor. These subsystems 
are connected to a local area network (LAN) along with a monitor and control (M&C) subsystem. The 
subsystems communicate with the M&C subsystem via a set of communications services implemented 
using manufacturing message specification (MMS) protocols. The subsystems communicate with one 
another by setting, sending, and reading MMS-variable structures. 

The TP- 13 receiver was integrated into the M&C LAN to support the Voyager telemetry downlink 
demonstration. Since most of the recent tracking activities at DSS 13 have not involved spacecraft, the 
telemetry processor was not at the station. Integrating TP-13 into DSS 13 involved writing an MMS- 
based interface for TP-13. At the same time, the TP-13 software was ported to the OS/2 operating 
system, requiring reimplementation of much of the low-level communication between the TP-13 software 
and hardware. 
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Fig. 3. DSS 13 configuration. 


The station monitor and control (SMC) subsystem provides the human operator with an interface 
to the subsystems shown in Fig. 3. As the name suggests, the SMC subsystem is used by the operator 
to monitor the status of the subsystems and to exercise control over them during a track activity. The 
LMCOA provides the capability to automatically execute operations procedures encoded in TDNs. In the 
SOE-driven automation demonstration, the DPLAN automatically generated the operations procedure 
in TDN format, and the LMCOA executed it. 


IV. Goal 1 : Automated Procedure Generation 

As part of the Voyager track demonstration, the capability to generate an executable operations 
procedure from high-level track information was demonstrated. The DPLAN planner uses high-level track 
information to determine appropriate steps, order constraints on these Steps, and determine parameters 
of the steps to achieve the high-level track goals given the equipment allocation. In generating the TDN, 
the planner uses information from several sources: 


(1) Project SOE — The project SOE specifies events from the mission/project perspective. 
As a result, the project SOE contains a great deal of information regarding the spacecraft 
state that is relevant to the DSN track, as well as a large amount of spacecraft information 
unrelated to DSN operations. Relevant information specified in the project SOE includes 
such items as the one-way light time (OWLT) to the spacecraft, notifications of the 
beginn in g and ending times of tracks, spacecraft data transmission bit-rate changes, 
modulation index changes, and carrier and subcarrier frequency changes. An excerpt 
from a Voyager SOE is shown in Fig. 4. 

(2) Project profile — This file specifies project-specific information regarding frequencies and 
pass types. For example, the project SOE might specify a high frequency, and the project 
profile would specify the exact frequency used. The project profile might also specify 
other signal parameters and default track types. An excerpt from a Voyager project 
profile file is presented in Fig. 5. 
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(3) LMCOA SOE (LSOE) — Using the project profile information, the project SOE is parsed 
into a more manageable form specifying only the information directly relevant to the 
specific pass. A sample LSOE for a Voyager telemetry downlink pass is shown in Fig. 6. 

(4) TDN knowledge base (KB) — The TDN knowledge base stores information on the TDN 
blocks available for DPLAN and LMCOA to use. This knowledge base includes informal 
tion regarding preconditions, postconditions, directives, and other aspects of the TDN 
blocks. It also includes information on how to expand the block parameters and name 
into the actual flat file entry in a TDN. 

(5) Equipment configuration — This details the types of equipment available and unique iden- 
tifiers to be used to specify the exact pieces of equipment to be used in the track. These 
include the antenna, antenna controller, the receiver, etc. 


*ID=DSNISOE52B1 S/C=032 BEG=94 360 094000 END=95 001 221500 
1032 42 00022 360 201500 ACE NCT S/C ANT STATUS, RCP/00.100/Y, HGA, RCP/00.100/Y, HGA 

1032 42 00022 360 201500 ACE NCT S/CTXR STATUS, OFF, EXC1/TWT2/1 8, 

LOW,EXC1/TWT2/18 
1032 42 00022 360 201500 ACE NCT 
1032 42 00022 360 201500 ACE NCT 
1032 42 00022 360 201500 ACE NCT 
!032 42 00022 360 201500 D42 BOT 

1032 42 00022 360 201 500 D42 ACD D/L, 1 -WAY, X 

1032 42 00031 361 095000 D42 EOT 


S/C RCVR STATUS, 2/S, ON 

S/C RNG MODE, OFF, OFF, ON, OFF 
S/C TLM X- MOD, 1 60CD, LOW, 61 DG 


Fig. 4. Sample Voyager project SOE. 


# switch IF3 to S-band (5), IF4 to S-band (7) 

IF3S 

IF4 X 

# LOW and HIGH subcarrier frequencies (kHz), from 870-14 
SUB_LOW 22500 

SUB_HIGH 360 
PRED_POWER 25.6 
SKY_FREQ 300000000 
CARRIER_FREQ 300000000 
CARRIER_UPDATE_RATE 500 
CARRIER_BANDWIDTH 1 
CARRIERJ.OOP II 
PC_NO 87.27 
SUB_UPDATE_RATE 100 
SUB_BANDWIDTH 1 
SUB_LOOP II 
PD_NO 275.42 

SYMBOL_UPDATE_RATE 100 
SYMBOL_BANDWIDTH 1 
SYMBOL_LOOP II 
ESJMO 0.871 


Fig. 5. Excerpt from Voyager project profile file. 



SKY_FREQ=300000000; 


PRED_POWER=25.6; 


95/005/16:30:00 CONFIG: IF3=X_KA; POLAR=RCP; 

PO_FILE=“/u/cm/lmcoa/lmcoa/src/kb/xv2_05.po";. 

95/005/16:30:00 TPR_PRECAL: BAND=X;. 

95/005/16:30:00 SOURCE: SRC=VGR2;. 

95/005/16:30:00 CARRIER: PC„NO=87.27; FREQ=300000000; BANDWIDTH=1.0; UPDATE_RATE=500; 
LOOP=2;. 

95/005/16:30:00 SUBCARRIER: PD_NO=275.42; FREQ=22500; BANDWIDTH=1.0; UPDATE_RATE=1 00; 
LOOP=2;. 

95/005/1 6:30:00 SYMBOL: ES_NO=0.871 ; BANDWIDTH=1 .0; UPDATE_RATE=1 00; LOOP=2;. 

95/005/16:30:00 ACQ_D/L: MODE=1-WAY; BAND=X;. 

95/005/16:30:00 TELEM: BIT_RATE=160; CODE_FLAG=1 ; MODJNDEX=61;. 

95/005/17:15:00 TELEM: BIT_RATE=2800; CODE_FLAG=1 
95/005/17:15:00 TELEM: MODJNDEX=72;. 

95/005/18:30:00 TELEM: BIT_RATE=160; CODE_FLAG=1;. 

95/005/18:30:00 TELEM: MODJNDEX=61;. 

95/005/20:00:00 END_OF_TRACK. 

Fig. 6. Voyager downlink telemetry pass LSOE. 



DSS 


Fig. 7. SOE-driven automation architecture. 

The overall flow of information in our approach to SOE-driven automation is shown in Fig. 7. In this 
approach, the project SOE is parsed into an LSOE that contains only that information relevant to the 
track of interest. DPLAN uses the LSOE and the project profile to determine the parameterized goals of 
the track. A sample set of track goals is shown in Fig. 8. This set of goals and a set of contextual facts 
derived from the LSOE are the planning inputs. The contextual facts are basically the planning system’s 
representation of the LSOE information. 

DPLAN reduces the high-level track goals into executable steps by applying knowledge about how 
to achieve specific combinations of track goals in the context of specific equipment combinations. This 
information is represented in the form of task reduction rules, which detail how a set of high-level goals 
can be reduced into a set of lower-level goals in a particular problem-solving context. Each task re- 
duction rule rigorously details the scope of its expertise in terms of track and equipment combinations. 
This information on the scope of applicability of the rule can be considered in terms of a track goal 
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(PO_required trackl) 
(spacecraft_mode_ohanges trackl) 
(track_goal spacecraft_track telemetry trackl) 
(track_goal spacecraft_track downlink trackl) 
(track_goal decode_data) 

(station-used trackl DSS13) 

Fig. 8. Sample tracking goals. 


hierarchy and equipment goal hierarchy, where the rule applies to all contexts below the rule in the 
relevant hierarchy (all specializations of its scope). 

Figure 9 shows a partial tracking-goal hierarchy involving the goals for telemetry, commanding, and 
ranging. Figure 10 shows a partial track hierarchy for antennas. A rule might specify how to achieve 
a tracking goal for a particular type of antenna. For example, a rule might specify how to achieve the 
telemetry tracking goal for a 34-m BWG antenna. Alternatively, a rule might specify a constraint on 
how to achieve telemetry when using high-efficiency antennas (HEFs). For example, a rule might specify 
that two receiver calibration steps, A and B, might have to be ordered so that A is always before B. 
By representing the track, equipment, and other hierarchies, the scope of various pieces of knowledge 
regarding antenna track activities can be naturally and easily represented. 

Using this problem specification, DPLAN then uses task reduction planning techniques, also called a 
hierarchical task network (HTN) [2,9], and operator-based planning techniques [11] to produce a param- 
eterized track-specific TDN to be used to conduct the track. The actual planner used to generate the 
TDN is a modified version of the task reduction planning component of the multimission Video Image 
Communication and Retrieval (VICAR) planner system [1] . 6 This track-specific TDN and the LSOE can 
then be used by LMCOA to operate the actual antenna to conduct the requested antenna track. 

For example, the task reduction rule shown in Fig. 11 states that in the equipment context of DSS 13 
and the tracking goal context of a downlink track, the abstract task of precalibration can be achieved 
by performing the lower-level tasks of inspecting the subsystems, connecting the subsystems, configuring 
the TPR, loading antenna predicts, and configuring the receiver. Furthermore, the subsystems must 
be inspected before connecting the subsystems, and so on. Also, some of these tasks are composite 
tasks and will be expanded later into more detailed tasks. For example, configuring the TPR involves 
configuring the IF switch, configuring the microwave switch controller (UWC)/TPR for precalibration, 
and performing the actual TPR precalibration. 


TELEMETRY COMMANDING RANGING 


TELEMETRY, COMMANDING TELEMETRY, RANGING 


TELEMETRY, COMMANDING, RANGING 
Fig. 9. Partial track hierarchy. 





6 S. Chien, R. Hill, and K. Fayyad, “Why Real-World Planning Is Difficult,” working notes of the 1994 Fall Symposium on 
Planning and Learning: On to Real Applications, New Orleans, Louisiana, November 1994. 
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EQUIPMENT SPECIFICATION 


GOAL TO ACHIEVE 


TRACKING GOAL SPECIFICATION 



Fig. 1 1 . Sample task reduction. 

Next consider the task reduction rule shown in Fig. 12. It states that in the equipment context of DSS 
13 and in any tracking goal context, the abstract goal of including a programmable oscillator (PO) can be 
achieved by adding the “load PO files” step followed by the “configure Doppler tuner” step. Additionally, 
these steps must be ordered with respect to the “connect subsystems” and “load antenna predicts” steps 
as indicated. 

In the context of specific tracks, generic configuration blocks will be specialized as appropriate to the 
details of the track. For example, in the equipment context of DSS 13, if the track context is downlink, 
telemetry, with symbol decoding requested, the task reduction rule in Fig. 13 states the specific receiver 
configuration block to use to configure the receiver. 

Considerable effort in computing the final TDN is devoted to determining the correct parameters 
(params) for blocks in the TDN. For example, Table 1 shows a configuration table used to determine the 
IF switch parameter for the TPR precalibration step. Depending on the communication bands used in the 
track, differing bands will be assigned to each of the communication pathways in the UWC/TPR. Based 
on the bands being used, the TPR IF switch parameter is also determined. This parameter setting is also 
determined during the decomposition process, so a correctly parameterized TDN can be constructed. 

The application of decomposition rules continues until all of the activity blocks in the TDN are 
grounded-that is to say, known blocks in the TDN KB can be used to instantiate the activities. This 
fully instantiated TDN can then be used with the LSOE by LMCOA to perform the actual track. Shown 
in Fig. 14 is the TDN for a Voyager downlink telemetry track using the PO for DSS 13. 


EQUIPMENT SPECIFICATION GOAL TO ACHIEVE 


CONNECT 

SUBSYSTEMS 



Fig. 12. Augmentation of the TDN required by an additional PO track goal. 


EQUIPMENT SPECIFICATION GOAL TO ACHIEVE 


TRACKING GOAL SPECIFICATION 


DSS 13 



DOWNLINK TELEMETRY DECODING 
TRACK TRACK SYMBOLS 


Fig. 13. A specialized task reduction. 
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Table 1 . TPR IF switch parameter determination. 


FYequency 

IF1 

IF2 

Parameter 

S-band 

* 

— 

1 

X-band 

* 

— 

1 

Ka-band 

— 

* 

2 

Q-band (45 GHz) 

— 

* 

2 

S- and X-band 

* 

* 

3 

X- and Ka-band 

* 

* 

3 




= COMMAND BLOCK = DECISION BLOCK 


= MANUAL TASK 


= TIME-TAGGED BLOCK 


Fig. 14. The TDN for a Voyager telemetry downlink track using the DSS-13 antenna 
(including telemetry mode changes). 


V. Goal 2: Automated Procedure Execution 

The previous section described how DPLAN automatically composes a procedure, called a TDN, 
from a knowledge base of procedure components, called TDN blocks, based on a request for a DSN 
service. In this section, we describe how a TDN is executed at a deep-space tracking station, DSS 13, 
using the LMCOA. As the name suggests, the purpose of the Link Monitor and Control Operator As- 
sistant is to assist DSN operators in performing tasks on the M&C subsystem. M&C tasks involve 
monitoring, configuring, operating, and troubleshooting deep-space communications equipment. The 
goal of automating the execution of procedures is to drastically reduce the amount of manual operations 
needed to support a DSN tracking activity by automating the execution of the operations procedures 
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for an activity. The benefits of achieving this, goal include reducing the operator’s workload, thereby 
allowing him or her to do other tasks, and reducing the amount of time required to perform a track. 


In order to understand better how LMCOA works, a short description of the operational context in 
the DSS-13 M&C subsystem is given first. Then we give an overview of LMCOA ’s architecture along 
with a description of how the TDN is executed. Finally, we discuss the results of the demonstration in 
which LMCOA executed the TDN for the Voyager telemetry downlink activity. 

A. Monitor and Control 

The M&C subsystem is the centralized site for monitoring and controlling the subsystems that comprise 
a deep-space communications link. The M&C console displays subsystem status data to the operator and 
also provides a graphical user interface for controlling the subsystems; it provides all the control functions 
needed for an operator to perform a tracking activity. For instance, Fig. 15 shows the M&C interface for 
performing the boresight procedure on the antenna subsystem. The state of the antenna is displayed as 
a table of values for azimuth, elevation, and offset. Control over the subsystem is exercised by pressing 
graphical user interface (GUI) buttons, such as the “MiniCal” button shown in the display. When the 
M&C subsystem is used without LMCOA to perform a tracking activity, the operator manually executes 
a sequence of actions (i.e., pushing buttons, typing data into data entry fields, switching displays, etc.) to 
satisfy the requirements of the activity. All control is exercised by the operator, who must also monitor 
the state of the subsystems in the communications link in order to close the control loop [5]. 
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Fig. 15. The M&C subsystem user interface for the boresight procedure. 
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B. Link Monitor and Control Operator Assistant 

There are several obvious drawbacks to operating the M&C subsystem manually. Certain DSN op- 
erations require continuous attendance by an operator over long periods of time, and some operations 
are highly repetitive and require large amounts of data entry. For instance, it is not unusual to conduct 
a KaAP track lasting 8 hours. During a KaAP track, the procedure called “load sidereal predicts” is 
repeated many times (see the KaAP TDN in Fig. 16 for an end-to-end view of the track); the “load 
sidereal predicts” procedure requires inputs by the operator each time it is conducted. We estimate that 
an 8-hour KaAP track requires about 900 operator inputs overall, if the track is performed manually. For 
this same track, under nominal conditions, LMCOA reduces the number of operator inputs to less than 
10. Since people tend to be more error prone than computers on simple repetitive tasks, it makes sense 
to assign these tasks to LMCOA, freeing the operator for the task of handling exceptions, which requires 
more intelligence and knowledge. 

The LMCOA performs the operations procedures for a tracking activity by executing a TDN, which is 
a procedure that is automatically generated by DPLAN, as described in the last section. DPLAN com- 
poses the TDN so that it contains the procedures (TDN blocks) needed for a specific tracking activity, 
and it orders them according to its knowledge of the dependencies that are defined among the blocks as well 



n 

= COMMAND BLOCK 

= MANUAL OPERATION 



= SUB-TDN 


o = DECISION BLOCK 


Fig. 16. A KaAP track TDN. 
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as by what it knows about the pre- and post-conditions of the blocks. The knowledge about interblock 
dependencies and about block pre- and post-conditions is passed to the LMCOA, whose task it is to 
execute the end-to-end procedure. The LMCOA receives the TDN in the form of a directed graph, where 
the precedence relations are specified by the nodes and arcs of the network. The blocks in the graph are 
partially ordered, meaning that some blocks may be executed in parallel. Temporal knowledge is also 
encoded in the TDN, which includes both absolute (e.g., “acquire the spacecraft at time 02:30:45”) and 
relative (e.g., “perform step y 5 min after step x”) temporal constraints. Conditional branches in the 
network are performed only under certain conditions. Optional paths are those that are not essential to 
the operation but may, for example, provide a higher level of confidence in the data resulting from a plan 
if performed. More details about TDNs are provided in [3]. 

To execute a TDN, LMCOA performs the following functions (see Fig. 17): 

(1) Load the parameterized TDN into the execution manager 

(2) Determine which TDN blocks are eligible for execution and spawn a process for each 
TDN block 

(3) Check whether the preconditions of each TDN block have been satisfied 

(4) Once the preconditions are satisfied, issue the TDN block commands 

(5) Verify that the commands had their intended effects on the equipment 

The operator interacts with the LMCOA by watching the execution of the TDN; the operator can 
pause or skip portions of the TDN, check the status of commands within individual blocks, and provide 
inputs where they are required. When part of a TDN fails, the LMCOA supports manual recovery by 
the operator by highlighting the point of failure and providing information about the preconditions or 
postconditions that failed to be satisfied. 


C. Demonstration Results 

Before LMCOA could perform a Voyager telemetry downlink activity, a number of steps had to be 
taken: 


(1) Knowledge was acquired about the operations procedures for performing this type of 
activity. This included collecting the commands for the activity; identifying each com- 
mand’s parameters and its source in the SOE and project profile files; identifying each 
command’s preconditions and postconditions and tying them to specific equipment state 
variables; making the commands into TDN block procedures; defining the precedence 
relations among the blocks; and displaying the status of the TDN’s execution to the 
operator. Though much of this knowledge had been previously captured for the subsys- 
tems at DSS 13, some TDN blocks had to be modified because the Voyager telemetry 
downlink activity is different from the previously performed activities. In addition, as 
described earlier in this article, a new subsystem, TP-13, was installed at DSS 13 for 
the demonstration. Thus, a significant amount of knowledge about TP-13 had to be 
acquired. 

(2) Portions of LMCOA had to be modified to accommodate the knowledge acquired for the 
Voyager telemetry downlink activity. In particular, the user interface had to be written 
for the new TDN, and rules were written in the situation manager. One of the long-term 
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Fig. 17. LMCOA architecture. 


168 





















goals for LMCOA is to generalise the architecture so that the TDN for any activity can 
be loaded and run without making modifications to any code, but for this demonstration, 
these modifications were necessary. 

(3) New MMS state variables had to be defined for the TP-13 subsystem so that it could 
communicate with the M&C subsystem and LMCOA. The MMS variable names are used 
by all of the subsystems on the M&C LAN, meaning that relatively minor modifications 
had to be made to TP-13, M&C, and LMCOA. 

(4) A TP-13 simulator was developed and installed in the network automation laboratory so 
that the communications interface could be tested between TP-13, M&C, and LMCOA. 
In addition, the TDN was debugged by executing it against the TP-13 simulator. The 
preconditions, postconditions, and time tags on TDN blocks often require a significant 
amount of testing before the TDN is validated. Such testing also helps identify missing 
knowledge. 

Once the preceding activities had been completed, the modified versions of LMCOA and TP-13 were 
installed at DSS 13 and the demonstration was conducted over the course of several weeks. 


VI. Knowledge Engineering 

The SOE-driven automation demonstration uses a knowledge base to store data about the domain and, 
as such, requires that knowledge engineering be performed in order to acquire and encode the knowledge. 
In this section, we describe the level and types of effort required to acquire, encode, and validate the 
knowledge base used in this demonstration. The levels of effort can be summarized as follows and are in 
units of 8-hour workdays. It took about 36 workdays to acquire the knowledge, almost the same amount 
of time to encode the knowledge, but three times as much time to validate the knowledge base. Table 2 
summarizes the knowledge engineering effort. Details about each phase of this effort are described below. 


Table 2. Knowledge engineering effort. 


Effort 

Time, 

days 

Acquire the knowledge 

36 

Background information 

14 

Reuse blocks from other TDNs 

4 

Interview experts 

11 

Produce TDN documentation 

7 

Encode the knowledge 

45 

Project SOE parser 

10 

M&C subsystem extended 

11 

Adapt LMCOA to Voyager TDN 

15 

Write, revise TDN flat file 

4 

Decomposition rules, syntax editor 

5 

Validate the knowledge base 

127 

TDN, pre- and post-conditions 

51 

M&C subsystem extensions 

19 

LMCOA adapted for Voyager TDN 

24 

Project SOE parser 

7 

Goldstone testing 

24 

DPLAN 

2 
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A. Acquiring the Knowledge 

The majority of knowledge represented in the system consists of the project SOE definition, the 
DPLAN knowledge base (used in generating the TDN), the TDN blocks themselves, the parameterization 
of the TDN blocks by both DPLAN and LMCOA, and the directives that make up the TDN blocks used in 
the demonstration. Much of this information comes from subsystem knowledge, for example, in defining 
the TDN blocks, identifying and understanding the subsystem directives within the blocks, and knowing 
how and when to parameterize subsystem directives within the blocks. Five subsystems were utilized for 
this TDN, two of which had never been used by LMCOA, and one that required more development and 
testing to incorporate it into the current TDN. Acquiring this knowledge took about 36 workdays. 

The information was obtained by several methods, including reviewing documents and learning about 
the existing software systems, interviewing experts familiar with a particular part of the domain, doc- 
umenting the knowledge, and participating in status meetings. Fourteen days were spent acquiring 
background information on DSN operations, SOEs, the M&C subsystem at DSS 13, and the LMCOA at 
DSS 13. 

A small but very valuable amount of time was spent reviewing the existing TDNs for DSS 13. A 
significant amount of knowledge was directly extracted from both the Ka-Band Link Experiment (KaBLE) 
and KaAP TDNs for use in the Voyager TDN. The operational unit of a TDN, namely a block, has proven 
advantageous in our previous work on TDNs. The reuse of blocks between different operations procedures 
is one key advantage. Twelve of the 22 blocks, or 55 percent, in the Voyager TDN came directly from 
the KaBLE and KaAP TDNs. Our plans for the next generation LMCOA include the use of a relational 
database to store TDNs, blocks, and their contents. A block that can be used in more than one TDN 
needs to be represented only once in the database. Changes could be made to one or more of these 
reusable blocks without having to make significant, if any, changes to each TDN requiring the block. 
Therefore, to answer the question put forth in the overview and objectives section, the effort to extend 
this demonstration to a different type of track can be greatly reduced, depending on the reuse of TDN 
blocks from existing TDNs. 

Most of the time, acquiring pre- and post-conditions occurred during the validation phase, due to the 
way the procedural knowledge is captured from operations personnel. Therefore, details about acquiring 
pre- and post-conditions are presented later in the validation section. 

Eleven days were spent interviewing experts in the areas of DSN operations in general, DSS-13 opera- 
tions for a spacecraft telemetry downlink pass, subsystem M&C, and project SOE definition. In addition, 
members of the Voyager project were interviewed in order to obtain parameters specific to Voyager that 
were required for parameterizing the TDN. Seven more days were spent generating graphical and written 
documentation of the TDN, and four days were spent participating in status meetings. 

B. Encoding the Knowledge 

The next step in the knowledge engineering process is to encode the knowledge into the right format 
in order for it to be machine readable. Due to the existing LMCOA, automated M&C subsystem, and 
planner, the majority of knowledge representations were designed and already used in previous efforts. 
In this effort, the majority of the encoding time was spent putting the acquired knowledge into a known 
format. Encoding the knowledge covers writing the TDN flat file, which specifies the Voyager TDN 
and is loaded by the LMCO; adapting the LMCOA to Voyager; writing the parser for the project SOE; 
extending the DSS-13 M&C subsystem to automate required subsystems; and developing DPLAN. Some 
highlights of the encoding effort are noted below. 

The DSS-13 M&C subsystem was expanded to include two new subsystems that were not required 
in other TDNs developed for the LMCOA. This included developing simple subsystem simulators for 
testing purposes. The integration of TP-13 into the M&C environment took the most amount of time. 
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The changes made to the LMCOA for the new Voyager TDN took 15 days. Some of this time was spent 
extending the definition of the LSOE to account for new items in the Voyager SOE. As other types of 
tracks and more subsystems are handled by the LMCOA, the definition of the LSOE will broaden to 
include operations (e.g., boresight, modulation index changes) required by different tracks. 

C. Validating the Knowledge Base 

This phase of the knowledge engineering required about three times the level of effort of either the 
acquisition or encoding phases. This phase includes on-site testing at Goldstone in the operational 
environment; testing and developing simulators in the DSN advanced systems laboratory; and validating 
the TDN, pre- and post-conditions, and other software modules. The validation phase took 127 workdays. 

Almost half of this time, 51 days, was spent validating the TDN and, especially, the pre- and post- 
conditions on TDN blocks. Validating the TDN involves making sure that the operations procedure 
is accurate. One reason this is a time-consuming task is that different experts have different ways of 
performing a part of the procedure or, over time, they revise and/or refine the operations procedure. 
This has been a constant in our experience developing TDNs. It is difficult to come to consensus on a 
single way to perform an operations procedure. After a TDN is in place and can be executed by the 
LMCOA, the operators can “see” the procedure more clearly and refine it. In any case, these reasons all 
point out the need to be able to modify TDN blocks and pre- and post-conditions. The format of the 
operations procedure must be simple so that it will be easy for developers and operations personnel to 
understand and maintain these procedures. 

Pre- and post-conditions were extremely time consuming to validate. Prom our previous experience 
building TDNs, we have found that operations personnel do not usually think in terms of pre- and post- 
conditions. The need for pre- and post-conditions is much more obvious after an initial demonstration 
of the baseline TDN. At this time, the operators observe subsystem actions occurring automatically in 
the TDN and identify when some actions occur before sufficient subsystem states have been reached. 
These states are then implemented as preconditions. Another reason why operators are not familiar 
with the concept of preconditions is related to how they operate the manual M&C environment. In this 
environment, the operators often have a lot of equipment preconfigured. A detailed question and answer 
session between the TDN developer and the operator proved useful for identifying what portions of the 
preconfiguration are preconditions for existing blocks in the TDN. Postconditions were more difficult for 
the operations personnel and developers to determine. 

The two main reasons for the large amount of time validating pre- and post-conditions are (1) the 
complexity of pre- and post-conditions and (2) subsystem support for them. Pre- and post-conditions 
can be very complex and, therefore, difficult to encode to have the desired effect. For example, a pre- or 
post-condition based on absolute or relative time can complicate the implementation of that condition. 
In an early version of the Voyager TDN, two TDN blocks occur in sequence, “start recording” and then 
“stop recording." The “start” recording block just tells the appropriate subsystem to start recording. 
The block then completes. The “stop recording” block does not execute until the appropriate time has 
been reached, according to the time in the SOE. Until this time, it appears that no blocks are executing 
in the LMCOA. That is true; however, the subsystems are recording data during this time. In order 
for the user interface to show that some activity was occurring, a postcondition was put on the “start 
recording” block. According to this postcondition, the “start recording” block will finish executing when 
the time has come to stop recording. Therefore, during the long data recording phase of the pass, the 
start recording block in the LMCOA remains in the executing state. 

Subsystem support is also required in order to make use of pre- and post-conditions. In this and 
previous TDNs, one or more subsystems did not provide status information that needed to be checked by 
a pre- or post-condition to the M&C subsystem. For example, in the Voyager TDN, files were downloaded 
to the Station Data Recorder (SDR). We had to modify this subsystem to send a “finished downloading” 
status back to M&C so that a postcondition could automatically test when downloading was completed. 
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Extensions to the M&C subsystem toot 19 days to complete. The majority of this time was spent 
incorporating TP-13, the receiver, into the M&C subsystem. New variables were added to the subsystem, 
and communications problems had to be overcome. This time also included writing subsystem simulators 
for testing M&C and the LMCOA in the DSN advanced systems laboratory. 

A significant amount of time, 24 days, was spent at DSS 13. Extensive testing and debugging that 
could not be performed in the DSN advanced systems laboratory environment were performed at Gold- 
stone. Much of this time was spent incorporating TP-13 into the DSS-13 M&C environment. Some of 
the difficulties encountered resulted in modifying sections of the knowledge base related to the M&C 
subsystem as well as the LMCOA. 

Finally, here are the levels of effort for debugging other software modules: Debugging the SOE parser 
took 7 days, and debugging DPLAN took 2 days. Another 24 days were required to debug the Voyager 
version of the LMCOA. 


VII. Discussion and Summary 

At the request of the TMOD Services Fulfillment Reengineering Team, a team from the network 
automation work area demonstrated SOE-driven automation for a Voyager telemetry downlink track 
at DSS 13 in February 1995. This SOE-driven automation demonstration accomplished the following 
significant results: 


(1) Demonstrated the ability to automatically generate procedures (TDNs) based on a 
project sequence of events (PSOE) using an artificial intelligence planning system called 
DPLAN 

(2) Demonstrated the capability for automated procedure (TDN) execution using the 
LMCOA system; this capability has been demonstrated before but never on a spacecraft- 
related track 

(3) Identified a set of gaps in the SOE definitions and project profiles that make it difficult 
to automatically generate operations procedures 

(4) Produced detailed measurements of the amount of effort required for knowledge engi- 
neering for necessary DPLAN and LMCOA knowledge bases 

(5) Showed that validating a TDN requires a lot of time and effort before it can be released 
for operational use 

(6) Pointed out the need for subsystems to report their status so that postcondition checking 
can be done more explicitly 

(7) Demonstrated that reusing TDNs can save a significant amount of effort in knowledge 
acquisition, knowledge engineering, and validation 
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